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Abstract 



This thesis investigates in the time domain a particular class of second order perturbations of a 
perfect fluid non-rotating compact star: those arising from the coupling between first order radial 
and non-radial perturbations. Radial perturbations of a non-rotating star, by themselves not emit- 
ting gravitational waves, produce a peculiar gravitational signal at non-linear order through the 
coupling with the non-radial perturbations. The information contained in this gravitational signal 
may be relevant for the interpretation of the astrophysical systems, e.g. proto-neutron stars and 
accreting matter on neutron stars, where both radial and non-radial oscillations are excited. Ex- 
pected non-linear effects in these systems are resonances, composition harmonics, energy trans- 
fers between various mode classes. 

The coupling problem has been treated by developing a gauge invariant formalism based on 
the 2-parameter perturbation theory (Sopuerta, Bruni and Gualtieri, 2004), where the radial and 
non-radial perturbations have been separately parameterized. Our approach is based on the gauge 
invariant formalism for non-radial perturbations on a time-dependent and spherically symmetric 
background introduced in Gerlach & Sengupta (1979) and Gundlach & M. Garcia (2000). It 
consists of further expanding the spherically symmetric and time-dependent spacetime in a static 
background and radial perturbations and working out the consequences of this expansion for the 
non-radial perturbations. As a result, the non-linear perturbations are described by quantities 
which are gauge invariant for second order gauge transformations where the radial gauge has 
been fixed. This method enables us to set up a boundary initial-value problem for studying the 
coupling between the radial pulsations and both the axial (Passamonti et al.,2006) and polar (Pas- 
samonti at el., 2004) non-radial oscillations. These non-linear perturbations obey inhomogeneous 
partial differential equations, where the structure of the differential operator is given by the pre- 
vious perturbative orders and the source terms are quadratic in the first order perturbations. In 
the exterior spacetime the sources vanish, thus the gravitational wave properties are completely 
described by the second order Zerilli and Regge- Wheeler functions. 

The dynamical and spectral properties of the non-linear oscillations have been studied with a 
numerical code based on finite differencing methods and standard explicit numerical algorithms. 
The main initial configuration we have considered is that of a first order differentially rotating and 
radially pulsating star, where the initial profile of the stationary axial velocity has been derived by 
expanding in tensor harmonics the relativistic j -constant rotation law. For this case we have found 
a new interesting gravitational signal, whose wave forms show a periodic signal which is driven 
by the radial pulsations through the sources. The spectra confirm this picture by showing that the 
radial normal modes are precisely mirrored in the gravitational signal at non-linear perturbative 
order. Moreover, a resonance effect is present when the frequencies of the radial pulsations are 
close to the first w-mode. For the stellar model considered in this thesis the gravitational waves 
related to the fourth radial overtone is about three orders of magnitude higher than that associated 
with the fundamental mode. We have also roughly estimated the damping times of the radial 
pulsations due to the non-linear gravitational emission. These values radically depend on the 
presence of resonances. For a 10 ms rotation period at the axis and 15 km differential parameter, 
the fundamental mode damps after about ten billion oscillation periods, while the fourth overtone 
after ten only. 
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Notations and convenctions 



• The signature of the metric is (—,+,+,+), thus time-like 4- vectors have negative norms. 

• The index notation of the tensor fields is the follows: Greek indices run from to 3, capital 
Latin indices from to 1, and small Latin indices from 3 to 4. 

• In a tensor expression we use the Einstein's sum convenction. 

• The Eulerian perturbation fields for the first and second order perturbations are denoted 
with the notation of the 2-parameter perturbation theory. Therefore, the upper index (1,0) 
denotes first order radial perturbations, (0,1) the first order non-radial perturbation and (1,1) 
their coupling. 

• In this thesis we have used the geometrical units in almost all the expressions. Thus, the 
speed of light and the gravitational constant are set G = c = 1. Therefore, we have that: 

1 s = 2.9979 x 10 5 km 
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Chapter 1 

Introduction 



Gravitational waves are the most elusive prediction of Einstein's theory of gravity. The indirect 
evidence of their existence relies on the observations of the binary pulsar PSR 1913+16 (Hulse 
and Taylor [55]), that shows a decay of the orbital period consistent with the loss of angular mo- 
mentum and energy due to the emission of gravitational waves. The prospect of starting a new 
astronomy based on gravitational radiation and providing a new corroboration of General Rela- 
tivity has motivated many theoretical and experimental researches. As a result, the detection of 
gravitational waves seems feasible in the next decade by an international network of Earth-based 
laser interferometer detectors (LIGO, VIRGO, TAMA300, and GEO600) 0, bar resonant an- 
tennas (EXPLORER, AURIGA, NAUTILUS, ALLEGRO) [Q and by the Laser Interferometer 
Space Antenna (LISA) [2[. Three scientific runs have been so far carried out by the LIGO detec- 
tors, in collaboration with GEO and TAMA detectors for two of the three runs and with the bar 
detector ALLEGRO for the last run. The data analysis of the first and second science run sets up- 
per limits on the gravitational signal emitted by a number of possible sources, such as stochastic 
background, coalescing binary stars, pulsars E]|31IIIl|71[n3IID- The third science runs have been 
performed with a higher sensitivity and the data analysis leads to a significant improvement of the 
gravitational radiation upper limits [9]. Meanwhile, a second generation of detectors is already 
in the design stage for the exploration of the high frequency band, up to several kHz (advanced 
GEO600 1971 . wide-band dual sphere detectors [21]), with an improvement of sensitivity up to 
two orders of magnitude with respect to the first-generation instruments. 

Gravitational radiation could provide new information about the nature of astrophysical sources 
and help in the interpretation of the dynamical evolution of many such systems. Among the many 
sources of gravitational waves, the oscillations of compact stars are considered of great interest 
by astrophysicists and nuclear physicists. The extreme conditions present in the core of compact 
stars make them a unique laboratory, where nearly all the modern theoretical areas of research in 
physics can be tested. 

In many astrophysical scenarios, compact stars may undergo oscillating phases. After violent 
events such as core collapse of a massive star, an accretion-induced-collapse of a white dwarf, 
or a binary white dwarf merger, the newly born protoneutron star is expected to pulsate non- 
linearly before various dissipative mechanisms damp the oscillations. Another system where 
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pulsating phases may occur is a massive meta-stable compact object, which is born after the 
merger of a binary neutron star system. The gravitational signal emitted by the stellar oscillations 
lies in the high frequency band (y > 1kHz) and strongly depends on the structure and physics 
of the star, for instance on the equation of state, rotation, crust, magnetic fields as well as on the 
presence of dissipative effects such as viscosity, shock formation, magnetic breaking, convective 
outer layers, etc. With a detailed analysis of the gravitational wave spectrum emitted by stellar 
pulsations we could infer through asteroseismology the fundamental parameters of neutron stars, 
such as mass, radius and rotation rate fi4l IT3"1 . This information is necessary for the nuclear 
physicists as they can test the equations of state proposed for the description of matter at supra- 
nuclear densities. However, the weakness of the gravitational signal and the noise associated 
with the location and technology of the detectors compels theorists to provide more and more 
accurate models to predict the spectral and wave form properties of the gravitational signal. These 
templates are indispensable for enhancing the chances of detection, by extracting the signal from 
the noise with statistical methods. 

The spectral and dynamical properties of the oscillations of compact stars have been exten- 
sively investigated during the last forty years in Newtonian and Einstein theories of gravity. Linear 
perturbative techniques are appropriate for the analysis of small amplitude pulsations both in the 
frequency fTTMIoT)!! and the time domain approach El EE EH EH E3- In General Relativity, the 
oscillation spectrum of a compact object, such as black holes and neutron stars, is characterized 
by a discrete set of quasi-normal modes (QNM). These modes have complex eigenfrequencies 
whose real part describes the oscillation frequency, and the imaginary part the damping time due 
to the emission of gravitational waves. The classification of QNM is well known for a large set of 
stellar models and can be divided schematically in fluid and spacetime modes. The fluid modes 
have a Newtonian counterpart and can be sub-classified by the nature of the restoring force that 
acts on the perturbed fluid element. The spacetime modes are purely relativistic and are due to 
the dynamical role assumed by the spacetime in General Relativity (more details are given in 
chapter |3] and reference therein). 

Rotating stars in General Relativity can be described with various approximations, such as the 
slow -rotation approximation ll52ll or recently with codes developed in numerical relativity 1371 . 
The former approach is based on a perturbative expansion of the equations in powers of the di- 
mensionless rotation parameter e = Q/Qk, where SI is the uniform angular rotation and is 
the Keplerian angular velocity, which is defined as the frequency of a particle in stable circular 
orbit at the circumference of a star. The measured period of the fastest rotating pulsar corresponds 
to a relatively small rotation parameter e ~ 0.3, which may suggest that the slow rotation approx- 
imation provides an accurate description of rotating stars even for high rotation rates. However, in 
these cases the accuracy of this perturbative approach is different for the various physical stellar 
quantities. For instance, the quadrupole moment shows an accuracy to better than twenty percent, 
while the radius of the corotating and counterrotating innermost stable circular orbits is accurate 
to better than one percent [18]. Nevertheless, a protoneutron star may be expected to have a 
higher rotation rate that is not possible to describe with the slow rotation approximation. These 
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regimes can be better addressed in numerical relativity by evolving the full set of non-linear Ein- 
stein equations l35lll05LfT0 3i. Furthermore, recent works on the core collapse ll33ll34ll . r-mode 
instability |94]|68]E], accretion from a companion [40], or supernova fall back material I113L 
show that a neutron star manifests a degree of differential rotation. 

These studies have clarified the effects of rotation on the dynamics of the oscillations, as well 
as showed the importance of using a relativistic treatment that takes into account the effects of 
the dragging of inertial frames. In particular, the presence in rotating stars of instabilities due to 
emission of gravitational waves has gained great attention. Almost all classes of oscillations of 
rotating stars, such as the f- and r-modes, are potentially unstable to the so-called Chandrasekhar- 
Friedman-Schutz (CFS) instability [24, 39]. This appears because beyond a critical value of the 
stellar angular velocity, a mode that in a corotating frame is retrograde and then has negative angu- 
lar momentum, may appear moving forward in the inertial frame of a distant observer. As a result, 
this inertial observer will detect gravitational waves with positive angular momentum emitted by 
this mode. Thus, the gravitational radiation removes angular momentum by the retrograde mode 
by making it increasingly negative and then leading to instabilities. The losses of the angular 
momentum through gravitational waves slow down the star on secular timescales; eventually the 
star rotates slower than a critical value and the mode becomes stable. These instabilities could be 
strong sources of gravitational radiation and also limit the rotation rate of neutron stars, providing 
a possible explanation for the measured rotation period of pulsars. Many studies are currently 
dedicated to understand whether viscosity, magnetic fields, shock waves on the stellar surface or 
non-linear dynamics of oscillations may saturate this instability. 

The non-linear analysis of stellar oscillations is more complex and only recently are some 
investigations being carried out, due to improvements achieved by the non-linear codes in nu- 
merical relativity [37]. Different methods have been used to investigate the properties of non- 
linear oscillations, such as for instance 3-dimensional general relativistic hydrodynamics code 
in Cowling or conformal flatness approximations [ l(3J[35]l, a combination of linear perturbative 
techniques with general relativistic hydrodynamics simulations [114], or a new method where 
the non-linear dynamics is studied as a deviation from a background, which is described by a 
stellar equilibrium configuration [ 102]. These works, which have been dedicated to investigate 
the non-linear dynamics of different astrophysical systems: non-linear oscillations of a torus or- 
biting a black hole 111411 . non-linear axisymmetric pulsations of uniform and differential rotating 
compact stars [35] and non-linear radial oscillations of non-rotating relativistic stars [ 1 02 1 , have 
revealed a new phenomenology associated with the non-linear regimes, the presence in the spec- 
trum of non-linear harmonics. These harmonics arise from the coupling between different classes 
of linear modes or from non-linear self-couplings [102], and have a characteristic that could be 
appealing for the detection of gravitational waves: their frequencies appear as linear combinations 
of the linear oscillation modes. Therefore, some of these non-linear harmonics (sub-harmonics) 
can emerge at lower frequencies than the related linear modes, and then be within the frequency 
range where the detectors have higher sensitivity. However, since the amplitude of non-linear 
perturbations is usually the product of the amplitudes of first order perturbations, in order to have 
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a detectable gravitational wave strain one needs non-linear effects that can enhance the gravita- 
tional signal, such as resonances, parameter amplifications or instabilities. 
Strong non-linear regimes are adequately studied with a fully non-perturbative approach. How- 
ever, many interesting physical effects of mild non-linear dynamics can be well addressed by sec- 
ond order perturbative techniques. An example is given by the analyses of black hole collisions in 
the so-called "close limit approximation", where the second order perturbations of Schwarzschild 
black holes have provided accurate results even in non-linear regimes where the pertur- 

bative methods are expected to fail. Non-linear perturbative methods have been successfully used 
also for studying linear perturbations of rotating stars, where the rotation is treated perturbatively 
with the slow rotation approximation I521I103I . 

An important aspect of second order perturbative analyses is that of providing an estimate of 
the error associated with the first order treatment, as there is not an a priori method to determine 
the accuracy of the linear perturbative results. Thus, the convergence and the corrections asso- 
ciated with any term of the perturbative series can be determined only by investigating higher 
perturbative orders. Furthermore, non-linear perturbative equations are usually a system of par- 
tial differential equations, thus their numerical integration is computationally less expensive than 
the full Einstein equations which are treated in numerical relativity. This relative simplicity of 
the perturbative approach may then provide accurate results and can also be used to test the full 
non-linear simulations. However, an extension to second order perturbative investigations is not 
always straightforward 1461 |4"T1 . Some issues may arise from the identification of the physical 
quantities among the second order perturbative fields or from the movement of the stellar surface 
in non-linear stellar oscillations ll01M102ll . 

In physical systems where the perturbative analysis can be described by more than a single pa- 
rameter, as for stellar oscillations of a slowly rotating star or mode coupling between linear pertur- 
bations, the multi-parameter relativistic perturbation theory fT9l flOO I can help the interpretation 
of the gauge issues of non-linear perturbations. The identification of gauge invariant quantities 
allows us to have direct information about the physical properties of the system under consid- 
eration. The construction of such quantities is not in general simple, but recent works l!53l 1541 
show how to build second order gauge-invariant perturbative fields from the knowledge of the 
associated first order gauge invariant perturbations. For a specific class of astrophysical systems 
a gauge invariant and coordinate independent formalism has been introduced nearly thirty years 
ago 1431 l431l for the analysis of one-parameter non-radial perturbations on a time dependent and 
spherically symmetric background. Recently, this formalism has been further developed [73 48], 
and has been used to study non-radial perturbations on a collapsing star 15H and for linear per- 
turbation on a static star [81]. 

The research project we have been working on aims to extend the perturbative analysis of 
compact stars at non-linear orders, in order to have a more comprehensive understanding of stel- 
lar oscillations and the related gravitational radiation. In particular, this thesis presents a gauge 
invariant formalism and a numerical code for studying the coupling between the radial and non- 
radial perturbations of a perfect fluid spherical star. The formalism for the polar perturbations has 
been worked out in a first paper 1591 . The formalism and applications to axial perturbations are 
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presented in [90]. Work in progress on the applications of the polar perturbative formalism will 
be presented in a future work. Radial and non-radial oscillations can be excited in the aftermath 
of a core collapse or by accreting matter on a neutron star. Radial perturbations of a non-rotating 
star are not damped by emission of gravitational radiation, but they can emit gravitational waves 
at non-linear order through the coupling with the non-radial perturbations. This picture changes 
in the presence of rotation, where the radial pulsations become sources of gravitational radiation 
and form a new class of modes, called quasi-radial modes. This coupling may be interesting for 
instance during a core bounce, where it is expected that an excitation prevalently of the quasi- 
radial and quadrupole modes. Even though the quadrupole component provides the dominant 
contribution to the gravitational radiation, the radial pulsations may store a considerable amount 
of kinetic energy and transfer a part of it to the non-radial perturbations. As a result, this non- 
linear interaction could produce a damping of the radial pulsations and an interesting gravitational 
signal. The strength of this signal depends naturally on the efficiency of the coupling, which is 
an effect worth exploring. In this thesis, we start to investigate this non-linear effect for small 
oscillations of a non-rotating star with the aim of including rotation in future works. 

The polar coupling, i.e. between the radial and the polar sector of non-radial perturbations, 
is expected to be a priori more effective than the axial case. Indeed, the linear polar modes have 
a richer spectrum than the axial sector and from the values of the frequencies and the damping 
times of the fluid QNMs, the resonances and composition harmonics should be more probable 
between the polar fluid and the radial modes. However, for the purposes of this thesis we have 
implemented a numerical code for the axial coupling for mainly two reasons: i) the axial coupling 
can have interesting physical effects, ii) the perturbative equations are simpler and enable us to 
understand better the issues related to the numerical stability and accuracy of the code as well as 
the effects of the low density regions near the stellar surface on the non-linear simulations, etc. 

When we consider the first order axial non-radial perturbations, we see that the only fluid 
perturbation is the axial velocity, which can be interpreted to describe a stationary differential 
rotation. Therefore, the linear axial gravitational signal does not have any dependence on the 
dynamics of the stellar matter. This picture changes at coupling order, where the differential ro- 
tation and first order metric perturbations can couple with the radial pulsations and source the 
axial gravitational waves. We will see in chapter |5] that this axial coupling produces a new class 
of quasi-radial modes, which can exist only for differentially rotating stars. 

This thesis is organized in seven chapters. In chapter EJ we introduce the perturbative for- 
malism used in this work, i.e. the multi-parameter relativistic perturbation theory and the gauge 
invariant formalism introduced by Gerlach and Sengupta and further developed by Gundlach and 
Martin Garcia (GSGM). In chapter [3] we describe the linear perturbations of a spherical star, i.e. 
the radial pulsations, the polar and axial non-radial perturbations. The equations for describing 
the coupling between the radial and axial and polar non-radial oscillations are presented in chap- 
ter |4j where in addition we discuss also the boundary conditions. In chapter |5j we present the 
proof of the gauge invariance of the perturbative tensor fields that describe the non-linear pertur- 
bations for this coupling. Chapter|6]is dedicated to the numerical code that simulates in the time 
domain the evolution of the coupling between the radial and axial non-radial perturbations. In 
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this chapter we give all the technical details relating to the code and the results of the simulations. 
Finally in chapter the conclusions and possible future developments are discussed. 

The appendix has seven sections. We have reported the source terms of the equations derived 
by Gundlach and Martin Garcia in section |X] In section |Bj we write the full expressions of 
the sound wave equation for the fluid variable H, while in sections [C] and |D] we respectively 
present the source terms of the perturbative equations that describe the coupling between the 
radial pulsations and polar and axial non-radial perturbations. In addition, in section |E| we give 
the tensor harmonics, while some of the numerical methods used in the numerical code are given 
in sections IFI and iGl 



Chapter 2 



Non-Linear Relativistic Perturbation 
Theory 

Exact solutions of the equations of physics may be obtained for only a limited class of problems. 
This aspect is particularly present in General Relativity, where the complexity of the astrophysical 
systems and the non-linear Einstein field equations allows us to describe exactly only simplified 
and highly symmetric cases. Among various approximation techniques, perturbation methods are 
appropriate whenever the problem under consideration closely resembles one which is exactly 
solvable. It assumes that the difference from the exactly solvable configuration is small and that 
one may deviate from it in a gradual fashion. Deviations of the physical quantities from their 
exact solutions are referred to as perturbations. Analytically, this is expressed by requiring that 
the perturbation be a continuous function of a parameter, measuring the strength of the perturba- 
tion. Although perturbative techniques are more appropriate for small values of the perturbative 
parameter, sometimes they can give reliable results also for mildly non-linear regimes as shown 
for example in the analysis of black-hole collision [92]. Hence, in many cases the validity limit of 
perturbative methods cannot be determined a priori. A more accurate estimation can be reached 
by studying the convergence of the perturbative series, which then involves the analysis of the 
second or higher perturbative orders. 

The gauge issue arises in General Relativity as in any other theory based on a principle of 
general covariance. The perturbative description of a physical system is not unique due to the 
presence of unphysical degrees of freedom related to the gauge, i.e. to the system of coordinates 
chosen for the analysis. This ambiguity can be eliminated either by fixing a particular gauge 
or by constructing perturbative variables which are invariant for any gauge transformation. In 
the former case, the properties and symmetries of the physical systems can help us to decide an 
appropriate gauge. In the latter approach, the identification of the gauge invariant fields is the 
difficult task. 

In this section we review the perturbative framework we have used for investigating the cou- 
pling between the radial and non-radial perturbations. In section |2~71 we report the main results of 
the multi-parameter perturbation theory introduced by Bruni et al. Q^i] and Sopuerta et al. [100]. 
Section l2~2l is dedicated to the formalism introduced by Gerlach and Sengupta [43||45l, which 
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has been further developed by Gundlach and M. Garcia [48 1, while in section l2~3l we outline 
the perturbative structure of our work which is based on the 2-parameter expansion of a static 
background. 

2.1 Multi-parameter perturbation theory 

Perturbation theory assumes the existence of two spacetimes, namely the background and per- 
turbed spacetimes. The former is the spacetime described by an exact solution of the field equa- 
tions, while the latter is the physical spacetime that the perturbation theory attempts to describe 
through the perturbation fields. The main requirement is that the physical description of the per- 
turbed spacetime slightly deviates from that of the background solution. 

Let us call M. the spacetime manifold. The multi-parameter relativistic perturbation formal- 
ism assumes the existence of a smooth multi-parameter family of spacetime models which are 
diffeomorphic to the physical spacetime: 

M X=( M > T X)- (2-1) 

The quantity T x represents a set of analytic tensor fields which are defined on M. and that de- 
scribe the physical and geometrical properties of the physical spacetime. The N-parameter vector 
A G R N labels any diffeomorphic representation of the physical spacetime Mr, and controls the 
deviation from the background quantities due to the contribution of the various perturbative pa- 
rameters of the system under consideration. In this notation the background manifold M& is then 
identified with the spacetime model Mf, = Mj. Furthermore, the validity of the Einstein field 
equations is assumed on any manifold M^: 

E{T x )=0. (2.2) 

In a perturbative approach, the comparison of perturbed and background variables is crucial for 
determining the accuracy of the perturbative description. In a physical theory based on a principle 
of general covariance as in General Relativity, this procedure requires a more precise definition 
that takes into account the gauge issue. Let us for instance consider a relation commonly used in 
perturbation theory, 

T x (q) = T b (p) + 5T X (p) , where 5T X « T b . (2.3) 

Here, q and T x are respectively a point and a tensor field in the perturbed manifold M?, while the 
point p and the tensors 7J, and 5T X belong to the background Mf,. From equation (12.31) the ten- 
sor Tx (q) can be considered as a small deviation of the background value %. However, we can 
notice that the perturbed and the background tensors are applied at points belonging to different 
manifolds. Therefore, in order to have a well posed relation it is necessary to establish a corre- 
spondence between these two points p and q and consequently between the three tensor fields in 
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the expression d2,3t . This point identification between the various representations of the physical 
spacetime is provided by a N-parameter group of diffeomorphisms T> v = |<£>^ : A G M w j, 

if : M x R N -> A*t (2.4) 

(p, A) i ^ (p(p,X) = ip^{p) , 

where the identity element corresponds to the null vector A = 0, i.e. tp^ (p) = p. The choice 
of the identification map ip is completely arbitrary and in perturbation theory this arbitrariness 
is called "gauge freedom". It is worth noticing that the gauge issue in perturbation theory is in 
general independent on the gauge related to the background spacetime, which fixes the system 
of coordinates only on the background manifold M^. In order to have a correct description of a 
physical system the physical observables have to be isolated from the gauge degrees of freedom. 
This can be accomplished by fixing a particular gauge, where the variables assume the correct 
physical meaning, or alternatively by determining a set of gauge invariant quantities. The latter 
procedure can be more difficult to realize, but it provides directly the physical quantities of the 
system. 

2.1.1 Taylor expansion 

A perturbative solution of the Einstein equations is built as a correction of the background solu- 
tion. This property, expressed in equation (I2.3t . allows us to approximate the physical variables 
and the field equation by their Taylor series. However, in order to define correctly a relativis- 
tic multi-parameter Taylor expansion some concepts related to the properties of the N-parameter 
group of diffeomorphisms have to be specified. 

In general, a consistent perturbative scheme should not depend on the order followed for per- 
forming two or more perturbations. We can then consider an Abelian group of diffeomorphisms 
which is defined by equation (I2.4t and the following composition rule: 



tp x o = ip x+>: , for y\,X'£R N . (2.5) 

Therefore, we can decompose every diffeomorphism ip? as a product of N one-parameter diffeor- 
morphisms: 

f\ = ^(AiA-,0) ^(o,.,o,A N ) • ( 2 - 6 ) 

It is also evident from the previous property and the commutation rules that there are AH equiv- 
alent decompositions of the diffeomorphism tp^. In any point of the perturbed manifold Mr, a 
diffeomorphism cp^ defines a N-parameter flow p^(p). This flow is generated by a vector field 
( G R N that acts on the tangent space of Mr-. In an Abelian group the generators of two different 



flows commute 
defined: 



= and a N-dimensional basis with the N independent vectors can be 
\(k} , where & = (0, .., ( k , ..0) . (2.7) 

l ) k=l. ' 



N 
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This basis generates the N one-parameter groups of diffeomorphisms of equation (I2.6t and can 
be used to decompose the vector field ( in its N components 



N 



C = J2 & where Cfc = (0, Cfc, -0) , (2.8) 



k=\ 

and to define the Lie derivative of an arbitrary tensor field T 



dipt T 



A fc 



where A fc = (0, .., A fc , .., 0) . (2.9) 

A fc =o 



The operator tpt is the pull-back associated with the flow ip^ . 



The Taylor expansion of the pull-back ip*. around the parameter A = is then defined as follows: 



n>0 



The last equality in equation d2.10t is a formal definition that will be very useful later for carrying 



out calculations with the Baker-Campbell-Hausdorff (BCH) formula. The definition (I2.10t and 
the diffeomorphism decomposition expressed in equation (12.61) allows us to define the Taylor 
expansion associated with the diffeomorphism <p?\ 

I N n p \ I N \ 

«V- E m^^^ EVik an) 

rti..njv>0 \p=l ^ / \P=1 / 

2.1.2 Perturbations 

In a particular gauge, the exact perturbations of a generic tensor field T are defined as follows: 

A7j = ^r x -r ff , (2.i2) 

where the tensors ATf and w*Tr as well as the background tensor Th are defined on the back- 

A A A 



ground spacetime A/5. The definition (I2.12t indicates that the background is the fundamental 
spacetime where all the perturbative fields are transported by the N-parameter flows and then 



compared with the background fields. The definition (I2.12t can be rewritten by using the Taylor 



expansion (12.1 It in the following form 

A7 ?- E (II?t1^-^ < 2 - 13 ) 



N V u p 



rti,..,njv>0 \p=l 



where the vector n = (n\, .., njv) controls the perturbation order of a tensor field with respect to 
the N-parameter, 

1 JV p=l 
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and riT = Ylp=i n v * s tne tota l perturbation order. 
2.1.3 Gauge transformations 

The representation of a perturbative field in the background manifold depends in general on the 
gauge choice (p?. Let consider two generic gauges represented by the diffeomorphisms ip and ip, 
which are generated respectively by the vector fields (^d, .., ^Ctv) an d (j Ci> ■•■>^Cn)- A gauge 
transformation is then defined by a diffeomorphism 

$ X :M^M (2.15) 

that for a given A, connects the physical descriptions determined in the two gauges tp and tp as 
follows: 

$ := p^ 1 o = <p_ % o ^ . (2. 16) 

The family of all diffeomorphisms that relate two gauges does not form in general a group, 

$:MxR N -» M (2.17) 
(p, A) ^ $(p, A) = $ x (p) . 



Since the gauge transformation <I>^ is a diffeomorphism we can extend to it some definitions used 
in the previous sections for the identification maps. For instance, the pull-back <!>*. of a generic 
tensor field T induced by the gauge transformation ( 12. 161 ) can be defined by using the expression 



(12.111) in the following way: 

= exp ^2 ^p£^(p S j ex P (~ ^2 Ap£^ P ^T. (2-18) 

The gauge transformation <1>^ is generated by the vector field which can be in general ex- 
pressed in terms of the two generators of the gauge transformation ip and i/j. In [ 100], the authors 
derive the relations between these gauge generators as well as among the perturbation fields by 
using the Baker-Campbell-Hausdorff (BCH) formula, which for two linear operators A, B is so 
defined: 

e A e B = e f(A,B) j (219) 

where the functional f(A, B) is given by, 

Pi qi Pm qm 

/( • m -&~*~ E ESfeTOTE^i ' <2 - 20) 

PuQi 
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and where in the previous expression the following notation for the commutation operators has 
been used [XiX 2 X 3 ■ -X n ] = [• • • [[Xx, X 2 ] ,X 3 ] , ■ ■ -,X n ]. By using the BCH formula, the 



gauge transformation (I2.18t reduces to a single exponential operator, 



$*.T = exp < 

A 



"1v,«jv>0 \p=l P 



(2.21) 



where 1 = £^ is the identity operator. 

The gauge transformations at every perturbative order are then determined with the following 
procedure: 



i) by using the definition ( I2.16t . one defines a new relation between the pull-backs ip*-T and ip*.T: 

A A 



A'A A 



(2.22) 



it) Thus, one can expand equation ( 12.221 ) by using the expressions ( 12. lit and A2.21I ). In doing 
that, one can use the linearity of the functional f(A,B) on the operators A and B and their 
commutators, and also the linearity of the operators A, B on the respective Lie derivatives. Hi) At 
the end, the desired relations can be determined by comparing the terms order by order (for more 
details see EH]). 

In case of 2-parameter perturbations N = 2, which is the parameter number used in this 
thesis, the gauge transformations at first order are the well known relations: 



(2.23) 
(2.24) 



At second order, the perturbation fields in the two gauges are related as follows: 



'£(2,0) 



'£(1,0) 



S^T + 2£^ 1) 5^T + [£e, n ^ + £ 



'£(0,2) 



'£(o,D 



£(i.i) T+ £ t s(°.i)r+ £ t 5 {1 ' 0) T 



+ 



°^£(i,i) { °^£(i,o) ' °^£(o,i) } 



(2.25) 
(2.26) 

(2.27) 



where {-, ■} is the anticommutator. Gauge transformations for higher perturbative orders can be 
found in reference [19]. 



2.2 Gauge invariant perturbative formalism (GSGM) 

Linear perturbations on a spherically symmetric background can be well described by using the 
formalism of Gerlach and Sengupta l43l RBI . With a 2+2 decomposition of the spacetime, the 
authors set up a covariant formalism to study linear non-radial perturbations on a time dependent 
and spherically symmetric background. Gundlach and Martm-Garcfa have further developed this 
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formalism for a self-gravitating perfect fluid f73ll48ll74l . The authors have specified a general 
fluid frame on which all the tensor fields and perturbative equations can be decomposed. This 
approach leads to a set of scalar gauge invariant fields and equations, which are easily adaptable to 
any coordinate system of the background. Hereafter we refer to this formalism with the acronym 
GSGM. 



2.2.1 The time dependent background 

The background manifold is a warped product M 2 x S 2 of a two dimensional Lorentzian man- 
ifold M 2 and the 2-sphere S 2 . The metric can be written as the semidirect product of a general 
Lorentzian metric on M 2 , g A B, and the unit curvature metric on S 2 , that we call 7 a &: 






' a /3 In 2 

1 " r 2 j ab 



(2.28) 



With Greek letters we denote the components defined in the 4-manifold, whereas the capital 
and small latin letters describe respectively the tensors in the M 2 and S 2 sub-manifolds. The 
scalar function r = r(x A ) is defined in M 2 , and can be chosen as the invariantly defined radial 
(area) coordinate of spherically-symmetric spacetimes. Besides the covariant derivative in the 
four dimensional spacetime, defined as usual 

ffa/3; M = , (2.29) 

we can introduced in the two submanifolds two distinct covariant derivatives 

9AB\C = , 7ab:c = , (2.30) 

where the vertical bar corresponds to the covariant derivative of M 2 and the semicolon to that of 
the 2-sphere S 2 . Moreover, we can introduce the completely antisymmetric covariant unit tensors 
on M 2 and on S 2 , eab and t ab respectively, in such a way that they satisfy: 

e A B\c = eab:c = 0, e AC e BC = -g 1 l, e ac e bc = j b a . (2.31) 

The energy-momentum tensor in a spherically symmetric spacetime has the same block diagonal 
structure as the metric, 

T ap = diag ( t AB , r 2 Q(x c ) lab ) , (2.32) 

where Q(x c ) is a function defined on M 2 . In this thesis we have used a perfect-fluid description 
of the stellar matter, thus we specialize the GSGM formalism to this case. Therefore, we have for 

T a /3, 

T a f3 = (p + p)u a Uf3 + pg a/3 , (2.33) 
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where p and p are the mass-energy density and pressure, and u a is the fluid velocity. The back- 
ground fluid velocity in spherical symmetry has vanishing tangential components, 

u a = {u A ,0). (2.34) 

The velocity ua and the space-like vector 

B A 

ua = —cabu => uau = , (2.35) 

provide an orthonormal two dimensional basis {ua, «b} for the submanifold M 2 . 

The metric g AB and the completely antisymmetric tensor eab can be written in terms of these 

frame vectors as follows 

9AB = -u A u B + n A n B , e A B = n A u B - u A n B , (2.36) 

while the energy-momentum tensor assumes the following form 

t A B = Puaub + pn A n B , Q = p. (2.37) 

In any given coordinate system for M 2 , {x A } , one can define the following quantity: 

V A = -r\ A ■ (2.38) 
r 1 

Then, any covariant derivative on the spacetime can be written in terms of the covariant derivatives 
on M 2 and S 2 , plus some terms due to the warp factor r 2 , which can be written in terms of va- 
The frame derivatives of a generic scalar function / are defined by 

/ = u A f\ A , f = n A f\ A , (2.39) 
which obey the following commutative relation: 

If )'-(/' )■=/"/' -^/- (2-40) 

Furthermore, a set of background scalars are introduced in order to write the background and 
perturbative equations in a scalar form: 

n = lnp, U = u A v A , W = n A v A , H = u A A , v = n A A . (2.41) 



Finally, the Einstein field equation for the background spacetime 

R 

Ra/3 - -w9a/3 = 8irT a/3 



(2.42) 
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in the 2+2 split are given by the following equations: 

- 2 (v A \ B + v A v B ) + (ivc + 3v c v c - ^\ g AB = 8vrf AB , (2.43) 

vff + v c v° - R = 8vrQ, (2.44) 

where R = \R^ is the Gauss curvature of M 2 . The conservation of the energy-momentum tensor 

T a/3 . Q = (2.45) 

leads to the energy conservation equation and to the relativistic generalization of the Euler equa- 
tion: 

Q+(l + ^(ji + 2U) = , (2.46) 
c 2 s n' +(l + l)v = , (2.47) 



and c 2 s is the speed of sound defined on the isentropic fluid trajectories as follows: 

2 d P 



(2.48) 



2.2.2 Perturbations 



Linear perturbations of a spherically-symmetric background can be decomposed in scalar, vector 
and tensor spherical harmonics. The perturbative variables are then completely decoupled in a part 
depending on the angular coordinates of the 2-sphere and a part defined on the submanifold M 2 . 
This expansion is really helpful because the perturbative problem is reduced to a 2-dimensional 
problem, usually a time and spatial coordinate. The tensor harmonics are decomposed in two 
different classes of basis, the so-called polar (even) parity and axial (odd) parity tensor harmonics. 
These transform differently under a parity transformation, namely as (—1)' for the polar and as 
(— for the axial. 

The basis for scalar function is given by the scalar spherical harmonics Y lm , which are eigen- 
functions of the covariant Laplacian on the sphere: 

^bylrn = _/(/ + l} Y lm . (2.49) 

where the integers (/, m) are respectively the multipole index and the azimuthal number. For a 
given I, the azimuthal number can assume the following 21 + 1 values: 



-I, -l + l, ..O.J-1, I. 
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There is not any axial basis for the scalar case. A basis of vector spherical harmonics (defined for 
I > 1) is 



Y 
S. 



lm 
a 

lm 



Y. 



lm 



polar , 



eX m axial, 



(2.50) 
(2.51) 



A basis of tensor spherical harmonics (defined for / > 2) for the polar case is given by 



Y 



lm 



Y lm lab , 



rvlm 



v lm , + ^) ylm,„, 
Y :ab H n Y 7a6 



(2.52) 



and for the axial class by the following definition: 



aim aim ■ oim 



(2.53) 



The explicit form of the tensor harmonics are given in appendix IE1 

The perturbations of the covariant metric and energy-momentum tensors can be expanded in 
the polar basis as 



/ idm vim 



j^pol lm ylm 



\ 



8i 



pol 

a/3 



h pol lm ylm r 2 ( R lm ylm + Qlm ylm^ 

( 5t%Y lm 



(2.54) 



r.pol lm v lm 
0L A x a 



V st 



pol lm ylm 



(2.55) 



r 2 5t 3 lm 7a6 Y lm + 5t 2 lm Z l a f J 



and axial basis as 



s 9 ; 



a/3 



Xj-CLX 



Lax lm aim 



i i,ax lm aim u I aim _i_ alm\ 
\ A n K^a.b ' ^bia) 

( o 



fi^ax lm cjlm 



(2.56) 



(2.57) 



\ 5tf lm S l ™ 5t lm (S%% + Sf£) 



In a spherically symmetric background the axial and polar perturbations are dynamically inde- 
pendent, and for a given multipole index I their dynamics does not depend on the value of the 
azimuthal number m. For simplicity, we can then study the non-radial perturbations on a spher- 
ical star by only considering the axisymmetric case m = 0. This approximation is not valid 
for instance in a rotating configuration, where axisymmetric and non-axisymmetric perturbations 
have different spectral and dynamical features. 

The true degrees of freedom on metric and matter perturbations can be determined by a set 
of gauge-invariant variables. In one parameter perturbation theory, see section 12.1.31 the first 
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order perturbation 5T of a generic tensor field T is gauge-invariant if and only if the following 
condition is satisfied [106]: 

L(T h = , (2.58) 

where % is the value of T on the background spacetime and £ is an arbitrary vector field that gen- 
erates the gauge transformation (see section l2~1.3l and reference [ 106]). By combining separately 
the polar perturbation fields /iab, h p A ° l , K, G, 5t A B, &t>Ai $t 2 , > and the axial ones h°f, h, 
5t A x , St , it is possible to determined the following set of gauge-invariant variables ll43ll45l. where 
for clarity the harmonic indices (I, m) are neglected, 
polar 

kAB = h A B - (pa\b + Pb\a) , (2-59) 
k = K-2v A PA , (2.60) 
Tab = 

^3 

T A = 

where Ta is defined for / > 1 

PA = 

axial 



~- 5t A B ~ t A B\CP ~ tACP\ B ~ t BCP\ A ' ( 2 .61) 

= 5t 3 -p c (Q lc + 2Qv c ) + l -^^-QG, (2.62) 

2 

= &tA-t A cP C -^QG {A , (2.63) 

= 5t 2 -r 2 QG, (2.64) 

, T 2 is defined for I > 2, and 

~ \r 2 G\ A ■ (2-65) 



ftA = hf ~ h\ A + 2\iva 
La = 5t a* - Qh°A , 
L = 5t-Qh, 



(2.66) 
(2.67) 
(2.68) 



where k A and La are defined for I > 1, and L for I > 2. Therefore, any linear perturbation 
of the spherically-symmetric background ( 12.281 can be written as a linear combination of these 
gauge-invariant quantities, which are tensor fields defined on the submanifold M 2 . The definition 
of the gauge invariant quantities of matter is valid for any energy-momentum tensor and not only 
in the case of a perfect fluid. 

The perturbations of a perfect fluid are given by four polar and one axial quantities. The polar 
velocity perturbation can be written as follows: 



Su a 



Y^ m oi m Y} rn 



(2.69) 
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where a is defined for I > 1 . The axial velocity perturbation is instead given by 

5u a = (oJ lm S l a m ^ 



(2.70) 



The functions 7 /m , a lm , (3 lm are defined on M 2 and describe the rate of the radial, tangential 
polar and tangential axial motion respectively. Furthermore, the axial perturbation [5 lm is gauge 
invariant for an odd-parity gauge transformation (see section l5'.l.lt . 

The mass-energy density and pressure perturbations can be written in the following form 
(using the barotropic equation of state) 



lm „xrlm 



bp = u> lm pV 



bp = c 2 s 5p . 



(2.71) 



In terms of these quantities it is possible to define a gauge-invariant set of fluid perturbations: 



a 

7 



a — p ub 



7 — n 



P B U A \B + ^U a (p B \ A - p A \ B ) 



(2.72) 
(2.73) 
(2.74) 



where in these expressions and in the remaining part of this section we do not write explicitly the 
harmonic indices (l,m). 

The gauge-invariant tensors ( I2.61l )-( l2~64l ). ( I2.67t and 12.681 ) for a generic energy-momentum 
tensor can be written in terms of the perfect fluid gauge invariant perturbations as follows: 

polar 



Tab 



1 



7 ( u A n B + n A u B ) + - (kAcUB + UAksc) u 



(P + P) 

+up (u A u B + c 2 s n A n B ) + P^ab , 
a(p + p)u A , 
pk + c 2 pu , 
0. 



c 



(2.75) 
(2.76) 
(2.77) 
(2.78) 



axial 



L A = P{p + p) u A 
L = 0. 



(2.79) 
(2.80) 



2.2.3 Perturbative equations 

The dynamics of linear oscillations of a time dependent and spherically symmetric spacetime 
is described by two independent classes of oscillations: the axial and polar perturbations. The 
perturbative equations can be expressed in terms of the gauge invariant GSGM quantities. In 
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addition, when a decomposition with respect to the vector basis {u A , n A } of the spacetime M2 is 
adopted they assume a scalar form 1481 . Here, we report the main procedure; see |48| for details. 
Polar sector: 

The tensor Icab can be decomposed on the frame {u A , n A }: 

kAB = v(~ u AUb + n A n B ) + §(u A u B + n A n B ) + ifj(u A nB + n A u B ) , (2.81) 
where 77, <fi and ip are scalars. It is useful to consider the following new scalar variable 

X = 4>~k + v, (2-82) 

in place of (/). Then, combining Einstein equations with the energy-momentum equations we can 
obtain the following set of equations: for I > 2, 

77 = 0, (2.83) 

for I > 1, 

-X + X " + 2( M -W = S x , (2.84) 

-k + c 2 s k" - 2c 2 .Uip' = S k , (2.85) 

= S^, (2.86) 

16ir{p + p)a = ip' + C a , (2.87) 

-a = S a , (2.88) 

+ -^j V = S u , (2.89) 

+ + = S 7 . (2.90) 

And finally, for / > 0, 

8^(p + p) 7 = (A0' + C 7 , (2.91) 

= -jfe" + 2J7V' + C w , . (2.92) 

where the expressions Of *S*^} ^ipj ^oti ^a? $oji 

can be found in Appendix lAl 

Axial sector: 

The perturbed Einstein and hydrodynamics equations can be written as 



- (I - 1) (Z + 2) r- 3 * = -16^e AB L A | B (2.93) 
P-c 2 s (n + 2U)p = 0, (2.94) 
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where, following references l43l l45ll we have introduced the gauge-invariant odd-parity master 
function ^ as 

y = r 3 e AB (r- 2 k A ) lB . (2.95) 

Once ^ is obtained as a solution of the odd-parity master equation d2.93t . the metric components 
k A can be recovered by means of the relation 

(I - l) (/ + 2 ) k A = 16vr r 2 L A - e AB (r ) |B . (2.96) 

The solutions are determined by specifying the initial values of f3, and ^ on a Cauchy surface. 
The fluid conservation equation ( 12.941 can be solved independently from the odd-master equation, 
as it depends only on the fluid perturbation f3. Its solution then provides a constant value of axial 
velocity j3 along the integral curves of u A [48 ]. However, in the odd-master equation $2.931 the 
velocity perturbation (3 couples with the background quantities, which being time dependent can 
source the non-radial oscillations. 

2.3 Non-linear perturbative framework 

The radial and non-radial perturbations are the two fundamental families of stellar oscillations, 
which have different properties with respect to the gravitational physics. In this section, we are 
going to investigate the main characteristics of the non-linear perturbations and their equations by 
adopting a two parameter perturbative scheme which allows us to distinguish at any perturbative 
order these two perturbation classes. The two parameter perturbation theory, is the N = 2 subcase 
of the multi parameter theory reported in section |2~T1 The parameter A denotes the family of radial 
perturbations, namely the class of polar perturbations with vanishing harmonic index / = 0. On 
the other hand, the second parameter e labels the class of non-radial perturbations with I > 1. 
With this notation the metric and energy-momentum tensors can be expanded as follows: 

o -5 4- Ap (1 ' 0) +^ (0 ' 1) 4- ^ p (2 ' 0) 4- A^ P (lll) 4- £2 a (0 ' 2) +<l(\ n f k ) f?Q7^ 

where the integers (n, k) are such that n + k > 2. The background tensors have been denoted 
with a bar, while the indices denote the perturbations of order i in A and j in e. A similar 
expansion can be done for the other fluid perturbations, i.e., velocity, mass-energy density and 
pressure. The equations at any perturbative order can be determined with a standard procedure: 

i) one introduces the perturbative expansions (I2.97l >. d2.98l > for the metric, energy momentum 
tensors and those related to the other fluid variables into the Einstein and conservation equations, 

ii) Taylor expand these equations with respect to the two perturbative parameters A and e, and 
eventually Hi) select the terms of the equation which refer order by order to the same perturbative 
parameter X n e k , where n, k G N. Let's carry out the analysis focusing on the Einstein equations, 
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the conservation equations can be addressed with the same method. We can start writing the full 
Einstein equations: 

E[g,il> A ] = G[g]-T[g,1> A ] = 0, (2.99) 
where G is the Einstein tensor, and ift A (A=l, . . .) the various fluid variables. After having in- 



troduced the perturbative expressions (I2.97l > and (I2.98l > into Eq. < I2.99I >. we obtain the following 
expression: 



E[g,^ A ) 



E b [g,i, A ]+XE^ 
A 2 



g 



(1,0)^(1,0) 



+ e 



#(o,i) 



g (0,D^(0,l) 



+ C_ E (W) [g(2.0) f ^.0)^(1,0)^ j(l,0) 



+ L J g (0,2) ^ ^(0,2) ^ F ( ,l) g, F (0,l) j + 0(A « > e fe) = o : 



(2.100) 



where the tensors J^ 1 ' ) and F^- ' 1 ' denote the set of metric and fluid variables of the radial 
and non-radial perturbations respectively. The linear differential operators E^'i' in the previous 
expression are defined as follows: 



rE 



(2.101) 



A i =eJ=0 



and Ef, = E^°'°\ They act linearly on the perturbation of order (i + j), and in general non- 
linearly on the background quantities. 

An interesting aspect of the second order perturbative equations is the presence in the ex- 
pansion (12.1001) of products between linear perturbations, which have been already determined 
by solving the first order perturbations. Therefore, in the non-linear perturbative equations they 
behave as source terms. The equation of order Ae can be then written as follows: 



#(l,i) 



(1.D 



j(i,o) g, F (o,i) 



(2.102) 



and a similar structure is also present in the A 2 and e 2 perturbative equations. The iterative proce- 
dure of the perturbation techniques implies that the part of the differential operators (I2.101t that 
acts linearly on the perturbations AV, as for instance in equation ( I2.102I ). is equal at any pertur- 
bative order. However, when the perturbative fields have different dependence on the coordinates 
the resulting systems of perturbative equations are different. This is the case for the radial pertur- 
bations, which unlike the non-radial do not have any angular dependence. In order to have more 
clear this distinction between the perturbative equations of the radial and non-radial perturbations 
we redefine the first order differential operators as follows: 



#(i,o) 



#(0,1) 



J NR ■ 



(2.103) 



where R and NR stand for "radial" and "non-radial" respectively. 
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Equation (12. lQQt has to be satisfied for arbitrary and relatively small values of the two per- 
turbation parameters (A, e). Therefore, each term of the expansion has to vanish and provide an 
independent system of equations associated with its perturbation parameters. The equilibrium 
configuration in this thesis is a spherically symmetric and perfect fluid star. The background 
spacetime is then determined by equations: 



Eb[g,j> A ] =0, 

which represent the Tolman-Oppenheimer-Volkoff (TOV) equations (see, e.g., [77 1). 
At first order we have two independent systems of equations for the radial 



(2.104) 



and for the non-radial perturbations: 

Lnr 



(1,0W1,0) 



0. 



(0,1)^(0,1) 



0. 



(2.105) 



(2.106) 



Lr 




= s 


jd,0) 


® 


Lnr 




= s 


jd.O) 


® F^ 1 )" 


Lnr 




= s 


2^(0,1) 





The second order perturbative equations instead can be divided in three independent systems of 
equations: the second order radial perturbation, the coupling between the radial and non-radial 
perturbations and the second order non-radial perturbations which are respectively given by the 
following expressions: 



(2.107) 
(2.108) 
(2.109) 



As we discussed above, the differential part Lr of the non-linear radial perturbative equations 
in ( 12.1071 is the same as in the first order equations (12.1051 . while the differential part Lnr 
of linear non-radial perturbative equations (I2.106t appears at second order for non-radial and 
coupling perturbations, respectively in equations ( 12. 1091 and (12.1081 . Perturbative tensor fields 
on a spherically symmetric background can be expanded in tensor harmonics (see section I2.2.21 . 
As a result, the angular dependence of the perturbations is decoupled from the dependence on the 
two remaining coordinates, which generally describe the time and a radial coordinate. Therefore, 
any perturbative tensor field 5T can be written as follows: 



JT(t,r,^0) = ^5T /m (t,r)W /m (0, 



(2.110) 



1(71 



where the quantity H lm denotes the appropriate tensor harmonics associated with the nature of 
the perturbative fields, i.e. scalar, vector or tensor as well as even or odd parity perturbation. The 
tensor ST lm (t, r) is the harmonic component of this expansion related to the harmonic indices 
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(l,m), which is determined by projecting the perturbation ST on the related tensor harmonic 
Tt lm (6, 4>) through the internal product associated with the 2-sphere S 2 : 



ST 



I'm 



(ST, 



n 



Im i 



T : H lm dQ 



(2.111) 



s 2 



where T : Ti has the following definition for a 2-rank tensor field T: 



Tnjlm rr n_/lm„,ac„ bd 
: H = Tabri cd 7 7 . 



(2.112) 



and 7 a fc is the unit metric of the 2-sphere S 2 . 

When we introduce the tensor harmonic expansion into the linear perturbative equations d2.105t 
and (12. 1066 . they assume the following expressions for the radial perturbations: 



„(i.o) ^(1,0) 

600 ' ^A.OO 



H 



00 



(2.113) 



and for the non-radial: 



E r r (0,1) , (0,1) 
l N r g\ m ' ^XiL 



l,m 

/00 _ v"00 



H lm = 0. 



(2.114) 



Since the spherical harmonic Tt = Y uu is a constant, in this section we will consider its value 
implicitly contained in the harmonic component T^' ^ of the radial perturbations. The decompo- 
sition in tensor harmonics allows us to describe independently the various harmonic components 
(Z, m) of the linear non-radial perturbations, where each component obeys the perturbative equa- 
tion (12.1 14t related to its indices (/, m). As we will see later this property is not in general valid 
in a second perturbative analysis. 

Equations (I2.107t - d2.109t that describe the non-linear perturbations assume the following form 
in terms of a spherical harmonic expansion: 
non-linear radial perturbations 



„(2,0) ,(2,0) 
600 



,(1,0) ,(1,0) 
00 ^ u 00 



(2.115) 



coupling radial/non-radial perturbations 



E T \ (1,1) , (1,1) 
Lnr zL'^ajL 
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Im 
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Im 



n 



1 in 



l,m 



non-linear non-radial perturbations 



J NR 



l,m 



(2,0) , (2,0) ]njlra_S^ V <7 \ F { ° A) 6?) 
%lm ' V>A, lm\ n - 2^ 2-> ^ |_ I'm' 63 * l"m" 
l',m'l",m" 



yA'm'jA"m" 



(2.116) 



(2.117) 



The presence of the source terms in the non-radial perturbative equations d2.117t prevents us 
from completely decoupling the perturbative components with different harmonic indices (l,m). 
In fact, the quadratic terms in the source couple different spherical harmonics according to the 
familiar law for addition of angular momentum in quantum mechanics. For instance, a complete 
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analysis of the quadrupolar case (I = 2) must take into account the source terms provided by the 
coupling of the indices (l 1 , 1") = (2,2) as well as in principle the indices (I', I") = (200, 198) 
and so on. Therefore, the dynamics of a second order perturbation T^f^ depends in principle 
on an infinite series of source terms F^^J ® fM, , which have to be solved by the related first 
order perturbative equations. In a non-linear analysis it is then crucial to select in the sources 
the dominant terms which provide the main contributions to the non-linear dynamics. This selec- 
tion is a standard procedure which has been used for instance in the perturbative analysis of the 
oscillations of a slowly rotating star, where the rotation has been treated perturbatively with the 
Hartle-Thorne slow rotation approximation l52ll54ll . In general, when the aim is the description 
of the gravitational radiation emitted by a physical system, the coupling between the quadrupole 
terms and the other moments with I close to I = 2 are expected to give the dominant contributions. 
In second order perturbations of Schwarzschild black holes 1^2"! 1351 . which have been used by the 
authors also for describing the collision of two BHs, the coupling between the quadupolar terms 
provides results which show good accuracy with respect to the non-linear simulations carried out 
in numerical relativity. 

In this thesis, we investigate the coupling between the radial/non-radial perturbations Ae, 
which obey the perturbative equations of the form <I2. 1 1 6b - In this case the equation can be 
easily decoupled, as for any harmonic component (I, m) of equation d2.116t the source terms 
contribute only with the following indices (l',m') = (0,0) and (l",m") = {1,0)- The source 
terms are then determined by solving at first order two system of equations, one for the radial 
perturbations (12.1131) and the other for the (I, 0) component of non-radial perturbations (I2.114I) . 

2.3.1 Time and frequency domain analysis 

The investigations of stellar oscillations are carried out in the time and frequency domain. These 
two approaches provide complementary information about the dynamics and the spectral proper- 
ties of the stellar perturbations. 

In the frequency domain, the time dependence of the perturbative fields is separated by the 
spatial coordinate, by assuming an harmonic dependence of the oscillations. Therefore, the tensor 
fields ( 12. 110b can be written as 

6T(t,r,0,<f>) = ^<5T im (r)e^W im (M), (2.118) 

lm 

where LOi m is in general a complex frequency associated with the harmonic indices (l,m). The 
different action of the radial and non-radial perturbations on the quadrupole of a spherical star, is 
also reflected on the mathematical nature of the frequencies u)\ m . A radially oscillating phase of 
a perfect fluid spherically symmetric star can be described by a Sturm-Liouville problem, whose 
solutions provide a complete set of normal modes with frequecies LOi m , where ui m G R. On 
the other hand, a non-radially oscillating dynamics is also a source of gravitational radiation, 
which damps the stellar oscillations. As a result, the non-radial spectrum is described by a set 
of quasi normal modes (QNM), where ui m are complex quantities whose real part describes the 
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oscillating frequency, and the imaginary part the damping time due to the gravitational emission. 
The spectrum of normal or quasi normal modes for radial and non-radial perturbations can be 
determined by eigenvalue problems, which can be set up by introducing the expressions (12.1181) 
into the radial and non-radial perturbative equations (12.1 13t and (I2.114t . The numerical methods 
to derive these results are for instance reviewed in references fTTHoH . 

In the time domain, there is not any harmonic assumption on the time dependence. The 
perturbative equations < I2. 1 13b - <T2. 1 14l > and d2.115t -( T2.117l ) are then integrated in a 1+1 numerical 
code, where one dimension describes the time and the other the spatial coordinate. This approach 
provides information about the time evolution of the perturbative variables of oscillating phases. 
In particular in gravitational physics researches, this method gives the propeties of the wave forms 
of the gravitational signal. In addition, the QNMs which have been excited in a time evolution 
can be determined by means of a Fast Fourier Transformation (FFT) of the time profiles. 



Chapter 3 

Linear Perturbations of Compact Stars 



Neutron stars oscillations have been extensively investigated with linear perturbative techniques 
both in Newtonian and relativistic approaches. The classical analysis of the oscillating star spec- 
trum has revealed the presence of various classes of modes which have been organized in a de- 
tailed classification. Linear perturbations are classified in two fundamental classes: the radial and 
non-radial oscillations. This definition respectively discerns the perturbations that have or not an 
angular dependence. In a non-rotating and spherically symmetric stellar model the adiabatic ra- 
dial pulsations are not damped by any dissipative or emitting mechanism. The single degree of 
freedom of radial perturbations, which represents the radial movement of the fluid, can be de- 
scribed by a Sturm-Liouville problem. Therefore, the radial spectrum is formed by a discrete and 
complete set of normal modes which provides a basis for decomposing the time evolution of any 
radially oscillating quantity by Fourier transformation. On the other hand, in relativistic stars the 
non-radial oscillations modify the stellar quadrupole and are damped by gravitational emission. 
The non-radial spectrum is then described by quasi-normal modes (QNM), which have complex 
eigenfrequencies where the real part provides the oscillation frequency of the modes, while the 
imaginary part identifies the damping time of the oscillations. 

The features of pulsation spectra are closely related to the properties of the stellar model 
adopted. The interpretation of the relations between the gravitational radiation and the source 
properties is the subject of the Astereoseismology, which is already a prolific area of the elec- 
tromagnetic astrophysics that has revealed important aspects of the internal dynamics of the Sun 
and non-compact stars. The high densities and strong physical conditions present in a relativistic 
stars prevent us from studying the neutron stars properties directly in Earth's laboratories. As a 
result, various equations of state have been proposed for describing the matter at supranuclear 
densities. An analysis of the gravitational spectrum related to these sources can settle this uncer- 
tainty, by determining the neutron star masses and radii with an accuracy sufficient to constrain 
the parameters of the equations of state proposed fT2l[T4ll . 

The general relativistic treatment of the radial pulsations started in 1964 with the work of 
Chandrasekhar ll23ll . The aim of these first studies was the stability issue of the stellar equilibrium 
configuration under radial pulsations. Subsequently, the interest moved to the investigation in the 
frequency domain of the spectrum features for various stellar models that are described with 
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more realistic equations of state (see [ 59 ] and references therein). The time evolution of the radial 
perturbations has been addressed quite recently in [95] and [101], in Eulerian and Lagrangian 
gauges. Ruoff in [95] has explored the numerical stability of the radial and non-radial oscillations 
when a polytropic equations of state of a star is replaced by a more realistic equation of state. 
Sperhake in [ 101], approaches the non-linear time evolution of radial pulsations of a polytropic 
non-rotating star in Eulerian and Lagrangian gauges. 

Non-radial oscillations of compact stars have been originally studied with the Newtonian 
theory of gravitation I112II . In this context, the gravitational radiation is due exclusively to the 
oscillations of the fluid, the emission rate is determined by the quadrupole formula [ 107 88 ] and 
the damping time by the expression E/E [16], where E is the pulsation energy and here the dot 
denotes the time derivative. Damping times of typical pulsation modes are very low [ 75 ] due to the 
weak coupling between matter and gravitational waves. The two classes of non-radial oscillations 
are the polar and axial perturbations. The axial perturbations have a degenerate spectrum, which 
is removed when the stellar model contains rotation, magnetic fields or non-zero stresses 11751 . 
For a perfect fluid non-rotating star the axial perturbations can describe a continuous differential 
rotation of the stellar fluid without any oscillating character. On the other hand, according to the 
nature of the restoring force that governs their dynamics the polar perturbations are classified in 
pressure (p), gravity (g) and fundamental if) modes. They have the following properties: 

f-mode: the fundamental mode is nearly independent of the internal structure of relativistic 
and Newtonian stars. It is the only mode present in the simplest stellar model, i.e. a zero 
temperature non-rotating star whose density is uniform. There is a single /-mode for any 
harmonic index I, and for a cold NS its frequency depends on the average density of the 
star. The f-mode reaches its maximum in amplitude at the stellar surface and does not have 
any node in the associated eigenfunction. Typical values of the frequencies and damping 
times are in the range 1.5 — 3.5 kHz and 0.1 — 0.5 s, 170*1 l32ll . 

p-modes: they are associated with the acoustic waves that propagate inside the star, where 
the pressure gradients act as restoring forces. A polytropic perfect fluid star is the simplest 
stellar model which can sustain these modes. The oscillating frequencies are higher than 
the /-mode, as they are related to the travel time of the acoustic wave across the star. The 
p-modes form for any harmonic index / a countable infinite discrete set, where the first 
element p\ has a typical frequency 5 — 6 kHz, damping time of one or few seconds and 
one node in the associated eigenfunction. The frequencies, damping time and the node 
number increase directly with the order of the mode. 

g-modes: These modes arise from temperature and composition gradients present inside 
the star. Gravity is the restoring force that acts through buoyancy forces. Like the p- 
modes, the g-modes form for any harmonic index I a countable infinite discrete set, but 
their frequencies are lower than the /-mode frequency and are inversely proportional to the 
order of the mode. The g-mode frequencies range from zero to a few hundred Hz, and in 
a perfect fluid star, which is the model adopted in this work, they are all degenerate at zero 
frequency. The typical damping time has an order of mag nitude of 10 6 s. 
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For more details about the mode classification see the monographs dedicated to this subject 1 1121 

EE 

The first relativistic analyses of non-radial perturbations of non-rotating stars is due to Thorne 
and his collaborators in a series of papers 1109111 1 ON 1 071 fT08l that date back to 1967. In General 
Relativity, the spectral properties and damping times of the stellar oscillations can be directly 
determined with eigenvalue problems, which provide the stellar QNMs. Subsequent researches 
have been dedicated to have a more complete understanding of the stellar QNMs, by extending the 
analysis to more realistic stellar models. For oscillations associated with the dynamics of the mat- 
ter variables (fluid-modes), the relativistic analyses provided some small corrections to the mode 
frequencies and more correct values of the damping times than the Newtonian approach [60]. 
However, the spacetime in General Relativity is not a static and "absolutum medium" on which 
the gravitational wave propagates, but has its own dynamical degree of freedom. This property 
adds to the Newtonian picture a new class of oscillation modes, namely the gravitational wave 
modes (w-modes) [23] which are high frequency and strongly damped modes that couple 
very weakly with the stellar fluid. This latter characteristic implies that the axial and polar space- 
time modes have similar properties. These purely relativistic modes can be separated in three 
classes. 

Curvature modes: are the standard u>-modes, which are present in all relativistic stars. They 
are associated with the "curvature bowl" present inside the compact star. The typical first 
curvature mode has frequency 5 — 12 kHz and a damping rate of tenth of milliseconds. For 
higher order w-modes the frequency increases and the damping rate is shorter. 

Trapped modes: Some of the curvature modes for increasingly compact stars (R < 3 M) 
can be trapped inside the potential barrier, when the surface of the star is inside the peak 
of the gravitational potential. The were first determined by Chandrasekhar and Ferrari [ 25 1 
for axial stellar perturbations. They have frequencies between a few hundred Hz and a few 
kHz and are more slowly damped by gravitational radiation than the curvature modes (few 
tenths of milliseconds). The spectrum of trapped modes is finite, the number of modes 
depends on the depth of the potential well and then on the compactness of the star. The 
main issue related to this class of modes is whether such an ultra-compact star can exist in 
nature. 

Interface or wji modes: were determined by Leins et al. [66|. There is a finite number of 
u>//-modes for any multipole I, which have frequencies that vary from 2 to 15 kHz and very 
short life (less than tenth of milliseconds). The existence of this family of spacetime modes 
may be associated with the discontinuity at the surface of the star. 

More details about the numerical techniques and physical properties of the QNM can be found in 
the reviews loTfllHBIIoH . 

In many works, the stellar fluid oscillations have been determined by neglecting the quantities 
associated with the gravitational field, i.e. the Newtonian gravitational potential or the metric 
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tensor of the spacetime. This method, known as the "Cowling approximation" E71 . provides 
frequencies and damping times of the fluid modes with an error usually less than 10 percent. 

The investigations of relativistic stellar perturbations as an initial value problem has been ad- 
dressed quite recently for spherical non-rotating stars in the context of gravitational collapse by 
Seidel [98], and for static stars by Kind, Ehlers and Schmidt [58]. This latter work has determined 
the set of perturbative polar equations and the appropriate boundary conditions for having a well 
posed Cauchy problem, which determine a unique solution. Subsequently, Allen et al. [11] and 
Ruoff [96 . 95 1, the latter using the ADM formalism, explored in the time domain the dynamics of 
linear polar non-radial oscillations of a non-rotating star for polytropic and more realistic equa- 
tions of state. These works provided important information about the numerical issues related to 
the numerical integration of the perturbative equations, and about the initial configurations which 
are able to excite the fluid and spacetime modes. By using the GSGM formalism, Nagar et al [81] 
extended the time domain analyses for investigating the non-radial perturbations of non-rotating 
stars induced by external objects, like point particles and accretion of matter from tori. 

Linear perturbations have been studied also for rotating relativistic stars with perturbative 
techniques and full non-linear codes. In this thesis we will study the non-linear oscillations of 
non-rotating stars, therefore we do not address here this subject. The interested reader can find 
accurate and up to date information in the review by Stergioulas [ 103 1. 

This chapter is organized in five sections. The equilibrium configuration is described in sec- 
tion while the background quantities of the GSGM formalism in section EOl The first order 
radial perturbations are introduced in section 13.31 while the polar non-radial perturbations are 
described in section |3~31 and the axial in section |3~51 

3.1 Background 

The equilibrium configuration is a non-rotating spherically symmetric star that is described by a 
static metric in Schwarzschild coordinates: 

g a pdx a dxP = -e 2 ^dt 2 + e 2A ^dr 2 + r 2 {d6 2 + sin 2 8d0 2 ) , (3.1) 

where the functions <&(r) and A(r) are two unknown functions that must be determined by the 
Einstein field equations. The radial coordinate r identifies for constant t and r a 2-dimensional 
sphere of area Anr 2 . 

The stellar matter is described by a single component perfect fluid, where by definition vis- 
cosity, heat conduction and anisotropic stresses are absent. This model, though simplistic, is 
suitable for a first investigation of non-linear oscillations and for a correct interpretation of the 
results. More realistic descriptions should consider the presence of magnetic fields, viscosity, 
crust, details of the stellar structure, superfluid and different particles, etc. We intend to take into 
account these specific elements in future works, also in order to avoid possible numerical insta- 
bilities which can always arise when the stellar model becomes more complex. The perfect fluid 
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energy-momentum tensor is given by the following expression: 

T a (3 = (p + p)u a Uf3 +Pg al 3, (3.2) 

where p and p denote the mass-energy density and the pressure in the rest-frame of fluid, and 
u a is the covariant velocity of the static background. The velocity is a timelike vector, thus its 
components can be derived by the normalization condition u a u a = —1. The covariant velocity 
assumes the following form: 

u a = (-e*, 0,0,0) . (3.3) 
The metric variable A can be related to a new function M by the following definition, 

r ( ~ -2A(r) 



M(r) = -[l-e-^) (3.4) 

In the stellar exterior this function assumes the constant value M = M(R S ), which is the gravi- 
tational mass of the star and R s is the stellar radius. In the Newtonian limit the functions M(r) 
and $(r) describe respectively the gravitational mass and the gravitational potential of the star. 

The Einstein (I2.42t and the fluid conservation equations d2.45t form a system of ordinary 
differential equations first derived by Tolman [111] and Oppenheimer and Volkoff [ 87 ] (TOV) in 
1939: 

/ M\ r 

= {^ + ^)7-2M (3 " 5) 



A r = ( Airpr - — ) - — — (3.6) 



M\ r 
7 1 ) r-2M 

P,r = -(p + p)$,r, (3-7) 

M r = 4vrpr 2 . (3.8) 

The integration of the TOV equations require the specification of an equation of state for the stellar 
matter p = p(p, s). We consider a cold neutron star at zero temperature. This approximation is 
certainly accurate for old isolated neutron stars in absence of accretion. A few seconds after a core 
collapse the temperature of a newly born neutron star rapidly decreases, and the thermal energy 
becomes much lower than the Fermi energy of the degenerate neutron fluid. The Fermi energy for 
nuclear densities p N = 2 x 10 u gcm~ 3 is about E F ~ SOMeV = 3 x 10 U K and increases for 
the supranuclear densities of the neutron star core. As a result, the thermal degrees of freedom can 
be considered frozen out. In this thesis we will investigate also the effects of coupling between 
radial pulsations and differential rotation, which is present within the first seconds of a proto- 
neutron star life. Therefore, the thermal and dissipative effects due to convective zones and shock 
formations near the stellar surface should have been included into the physical model. In order 
not to complicate our investigation of the non-linear oscillations we neglect these effects with the 
aim to include them in future works. 
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In the present work, the star is then described by a barotropic fluid p = p(p), which is param- 
eterized by a polytropic equation of state (EOS): 

p = kp T , (3.9) 

where k is the adiabatic constant and T the adiabatic index. The background speed of sound in 
the fluid is then given by 

c 2 = dp/dp. (3.10) 

The TOV and fluid equation of state provide a one-parameter family of solutions that depend on 
the stellar central density p c . Its numerical integration is described in section lo"2l 

The exterior of a non-rotating star is a Schwarzschild spacetime represented by the following 
line-element in Schwarzschild coordinate: 



ds z = - ^1 - — J di A + ^1 - — J dr 2 + r 2 d6 2 + r 2 sin 2 6 d<j) 2 (3.11) 

where M is the gravitational mass of the star M = M(R S ). The internal and external solutions 
have to match on the stellar surface. Therefore, the following condition must be satisfied by the 
metric variable 3>, 

$(R S ) = -A(R S ) = ~ In (l - . (3.12) 

3.2 GSGM background quantities for linear perturbations 

The linear non-radial perturbations on a non-rotating star can be studied with the gauge-invariant 
perturbative formalism set up by Gerlach and Sengupta and further developed by Gundlach and 
Martm-Garcfa, which has been introduced in section 12.21 This formalism can be specialized to 
the case of a static background by choosing the static frame vector basis {u A ,n A }, 

u A = ( e -*,0), n A = (0,e- A ). (3.13) 



The associated background scalars d2.41t assumes the following expressions: 

e" A 

fl = U = 0, zy = e - A $ r , W = , (3.14) 

r 

while the frame derivatives of a generic scalar function f(x A ) are given by: 

f = e~*f,t f' = e~ A f, r (3.15) 



3.3 Radial perturbations 

The radial perturbations of non-rotating stars preserve the spherical symmetry of their equilibrium 
configuration. The geometrical and physical quantities of a star deviate from the background 
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values only along the radial coordinate. Therefore, the validity of Birkhoff 's theorem implies that 
for a non-rotating star an external observer cannot receive any gravitational information about the 
pulsating stellar dynamics. Adiabatic radial oscillations can be described in terms of their normal 
modes. This characteristic is more evident in the frequency domain analysis, where the radial 
perturbative equations can be set as a Sturm-Liouviile problem, which provides a complete and 
discrete set of the eigenfrequencies of the normal modes and their associated eigenfunctions. 

The time domain analysis of the radial oscillations has been investigated in various gauges 
and with different set of variables and equations 1771 l95l l48l . The most common perturbative 
quantity used to describe the unique degree of freedom of radial perturbations is the Lagrangian 
displacement ^(x a ). This is a vector field that provides at any instant of time the position of 
a fluid element with respect to its equilibrium position. The Lagrangian displacement obeys a 
wave equation that can be studied in the time or frequency domain. In this work, we describe the 
radial perturbations by using the GSGM formalism. This choice allows us to use a more uniform 
set of perturbative variables for both the radial and non-radial oscillations. The Lagrangian dis- 
placement and its eigenvalue equation will be useful later for setting up the initial configuration 
(16.3.31 and for estimating the movement of the surface along the evolution d6.3.4t . However, it 
is worth noticing that the GSGM formalism fails to be gauge invariant for radial perturbations, 
for instance some of the polar gauge-invariant tensors d2.59t -( IT63l cannot be even defined. This 
property does not produce any limitation to the approach of this thesis, as we will prove later that 
the non-linear perturbations that describe the coupling will be gauge-invariant only for a fixed 
radial gauge. 

The metric of radial perturbations is given by the following expression: 



/ \ 



x (1.0) 



(3.16) 



V r^K^hab ) 



where the constant value of the harmonic scalar functions Yqo is implicitly contained into the 
perturbative variables, while the tensor h A 's assumes the form: 

h AB = V {lfi h A B + 4> m (u A u B + n A n B ) + ^°\u A n B + n A u B ) . (3.17) 

The four scalar quantities rj^ 1 ' ' ,(f>^ 1 ' \t/j^' ' and K^ 1 ' ' are functions of the coordinates (t,r) 
and describe the metric radial perturbations. A gauge choice that considerably simplifies the 
perturbative equations is the following: 

^(i,o) =0> fc(i,o) =0j (3.18) 
which makes the radial metric diagonal, 



e 2* ( x (i,o) _ 2r? (i,o)) 

e 2A x (l,0) 





6g% 0) = e 2A X (1 ' 0) I . (3.19) 
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where we have used the definition ( 12.821) for the new metric variable x^ 1 ' ^' 

x (i,o) =77 (i,o) + 0(i,o) _ (3 _ 20) 

As shown in reference [74], the relations (13- 1 8I > do not fix completely the gauge of radial pertur- 
bations. There is a residual gauge degree of freedom that can be used in the boundary conditions 
to impose the vanishing of the metric perturbation r/ 1 ' ) on the stellar surface. 
The radial perturbation of the fluid velocity is given by equation (12.691 . 

^' 0) = ((4^- r?(1 '° ) ) ^ eV^Uo) , (3.21) 

where the scalar function 'y( 1 ' ) depends on the (t, r) coordinates and describes the fluid element 
velocity along the radial coordinate. The pressure and density perturbations are for a barotropic 
fluid given by: 

Sp M)=JW)p, 5p^=c 2 s Sp^°K (3.22) 
3.3.1 Radial perturbative equations 

The perturbations of Einstein and conservation equations lead to a set of four partially differential 
equations for the four variables ^v 1 ' ),^ 1 * ),^ 1 ' ),^ 1 ' ), see |4*81 . Before writing the system 
of equations, it is more convenient for calculation purposes to adopt a slightly different set of 
radial perturbations. The density perturbation u/ 1 ' ) will be replaced by the enthalpy perturbation 
H^'°\ while the metric quantity x^ 1,0 ^ W11 l be changed with the variable S^'°\ in order to use a 
set of variables consistent with the one we will use for the non-radial polar perturbations. These 
two new perturbations are defined as follows: 

H m ^8jW = %p_ w(1|0) gm ^ 

p+p p+P r 
The dynamics of the four radial perturbations S^- 1 ' ' , r/ 1 ' ) , 

i,o) )7 (i,o) is 

then governed by 

the following set of three partially differential equations that does not contain the quantity t^ 1 ' ', 



-4vr(p + p)r] e* +A 7 (1 ' 0) , (3.24) 
-e*" A -An(p + p)r e*+ A - (^npr 2 + ±) e * +A S^ , 

(3.25) 

S ( J' 0) = -8tt(p + p) e *+ A 7 ( 1 .°). (3.26) 



i J Uvrpr + -s- ] + -e 2A 
ci J \ r z J r 



(i-o) 
7t 
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Equation (I3.26t has been derived by using equation (34) in Ref. 1731 . The remaining metric 
variable t/ 1 - ) is then obtained by a constraint, which is the following elliptic equation: 



^(1.°) = 4ir(p+p)r 



.#(1,0) + f 1 + £Td,o) 



„2A 



(3.27) 



The solutions of radial pulsations must satisfy the Hamiltonian constraint on the initial Cauchy 
hypersurface and all along the evolution. The equation of this constraint is given by: 



Si 1 ' ) 



■2A 



Q _ 2 2M\ a(10) 
8npr 1 ^ ) 



8vr 



E±l H m 



(3.28) 



The system of equations (I3.24t - (l3.28t is equivalent to the one used by Ruoff [95], cf. also [77]. 
The evolution equations (13.241) and (13.251) can be combined in order to determine a wave equation 
for each of the variable involved flv 1 ' ) or The resulting equation for the enthalpy requires 

another equation to close the system, as it contains the metric perturbation S^ 1,0 ) as one of its 
terms. On the other hand, the radial velocity ^t 1 ' ) satisfies a single wave equation. Therefore, 
the radial velocity 7( 1,0 ) can be used to represent as well as the Lagrangian displacement the 
single degree of freedom present in a radially oscillating configuration. This equation can be 
determined by differentiating the two equations ( 13.241 ) and 03.25b with respect to the radial and 
time coordinate respectively, then we can take an appropriate linear combination of them and 
introduce some of the radial equations in it. This procedure leads to the following wave equation: 



7 J 0) - cl e 2 (*" A ) ^ + d x {r) 7 ^ 0) + d 2 (r) = , 



(3.29) 



where the background coefficients d\{r) and d 2 {r) are the following: 



-2$ 



d 2 (r) 



i , M \ 1 dc 2 s 

(„ + P )(4apr +pr )^4 



M 2 

Air (p - 2p) r + — - - 



M 

+ 47T£> r H J 



I , M 

[p + p) Airp r + ~2 



2 6M . . / ivj. \ OA 

IQ-irp - -2 + -^3 8irr {p + p) [ iirpr - — ) e 



M 
r 

M 



1 3d 



(3.30) 



2 4\ 



(3.31) 



Having solved equation ( 13-291 ) for we may use the first order evolution equations to deter- 

mine the enthalpy H^'°\ and the metric S^ 1 ' ) = ^( 1,0 )/r and j/ 1 ' ) variables. 



3.3.2 Boundary conditions for radial perturbations 

The physical solutions of the radial perturbation problem have to satisfy the boundary conditions 
at the origin and surface. The origin of coordinates r = must be a regular point for the pertur- 
bative fields and equations. The analysis of the Taylor expansion of the perturbative fields around 
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the centre leads to the following expressions: 

g (l,o) = S^°\t)r + 0(r 3 ), (3.32) 
^(i.o) = ^(1,0) ^ + 0(r , 2) ^ (3 _ 33) 

#(1,0) = H ^\t)+0{r 2 ), (3.34) 
7 (i,o) = 7 (i.o)(t) r + o( r 3). (3.35) 

The physical condition on the stellar surface is the vanishing of the total pressure on the perturbed 
surface: 

E = {x^eM | P (x fl )=0} (3.36) 

In a perturbative approach, the small movement of the surface from its equilibrium position can 
be described by the Lagrangian displacement £(i,o) = £ r 0& a ) that gives the position of the 
perturbed surface with respect to the background surface. Thus, the total pressure (I3.36t can be 
Taylor expanded around the equilibrium surface 

p (a£ + A e) = P (4) + ^P im (4.) + 0(X 2 ) , (3.37) 

where x^, are the coordinates of the background surface, A controls the strength of the radial 
perturbations and Ap^ 1 ' ^ is the radial Lagrangian perturbation of the pressure. This last quantity 
can be written in terms of the Eulerian perturbation with the well known relation that connects 
the perturbations in these two gauges (see e.g. Il771l ). 

Ap (l,o) =6p m + £ kl0) p. (3.38) 

The condition (13.361 and the perturbative expansion d3.37t leads to the vanishing of the back- 
ground pressure p (x^) = and its Lagrangian perturbation Ap^ 1 ' ' = on the equilibrium 
surface of the star. 

The condition p (x^) = will be imposed for the numerical integration of the TOV equa- 
tions. The Lagrangian perturbation of the pressure can be written in terms of the MTW radial 
renormalized displacement function Q as 

r 2 A (1,0) = _< p + p ) d 2 (3 39) 

or 



where ( is so defined: 



However, we cannot use directly the expression d3.40t . as the Lagrangian displacement is not 
a dynamical variable in our set of radial perturbations. The rate of the surface movement is 
described instead by the radial component of the fluid velocity, namely 7c 1 ' ) which is related to 
( as follows 

Ct = r 2 e- A 7 ( 1 '°) . (3.41) 
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By differentiating equation (I3.39t with respect to the time one finds the following boundary con- 
dition on the surface for ^t 1 ' ): 



(p + p)cle-*( 



r 2 e -A 7 (l,0) 



0. 



(3.42) 



In a polytropic equation of state the pressure, mass-energy density and the speed of sound vanishes 
on the static surface, thus the equation (I3.42t implies that the radial variable -y^ 1 ' ) and its spatial 
derivative 7^ must be finite. On the other hand, when the mass-energy density and the speed 
of sound are non-null on the surface the following condition on -yt 1 ' ) is required 



r 2 e -A 7 (l,0) 



0. 



(3.43) 



r=R 



This is the only physical surface condition. The behaviour of the metric perturbation S^ 1 ' ) and 
the enthalpy f/A 1 ' ) can be directly deduced by the perturbative equations d3.26t and (I3.24t respec- 
tively. For the other metric perturbation 7/ 1 ' ) one can use the residual gauge degree of freedom 
that has not been fixed by the radial gauge (13.181 . see reference 1741 . In accordance with the phys- 
ical properties of the radial perturbations, the more appropriate choice is a null value of ry( 1,0 ) on 
the surface. This allows us to eliminate all the gauge fields present in the external spacetime f74l . 



3.3.3 Frequency domain analysis of radial perturbations 



The time domain integration of radial perturbative equations (I3.24t - fl3.26t requires the choice of 
the initial values for the radial variables. A method used to determine the initial profile of one of 
the radial perturbations is to provide the eigenfunction associated with the particular radial mode. 
This method allow us to select and excite the radial modes of our interest and then simplify the in- 
terpretation of the non-linear effects due to the coupling with the non-radial oscillations. In order 
to determine the eigenfunctions for the radial variables ^i- ) we can elaborate the most common 
eigenvalue equation used in literature 11771 59 1, which is the wave equation of the Lagrangian 
fluid displacement ^' \ By using the renormalized Lagrangian displacement ^v 1 ' ) defined in 
equation (I3.40t the wave equation reads [ 77 1 : 



0, 



(3.44) 



where W, P, Q are functions of the radial coordinate r only. They are defined as: 



r 2 W 
r 2 P 



r 2 Q 



(p + P) e 3A+ *, 



(P + P) 



$ 2 



r 



8irpe 



2A 



Q A+3$ 



(3.45) 
(3.46) 

(3.47) 
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The radial component of the velocity perturbation ^yt 1 ' *) obeys the wave equation (I3.29t . An 
eigenvalue problem can be then set up by introducing the harmonic ansatz (12.1 18t in the vari- 
able However, in order to simplify the analysis the wave equation (I3.29t can be written 
after some manipulations with the same form of equation (I3.44l >: 



0, 



(3.48) 



where we find it convenient to use the variable y^'°\ which is related to c^ 1,0 ) by 

j,(l,0) =r 2 e -A 7 (l,0) =C (1.0) =r 2 e -* e 



(3.49) 



and W, P, Q are the same three functions defined in equations (I3.45l) - J3.47l) . We notice that 



the wave equation (I3.44t can be also determined by time differentiating the equation (13.441) and 



then using the relation (I3.41t that connects the two variables c^ 1 ' ) and 



7 



(1,0) 



However, the 



definition (13.411) shows that the solution of the two equations d3.44t and d3.41t disagree at most 
for a function that depends on the radial coordinate r only. This function can always set to zero 
with an appropriate choice of the initial conditions. 

The eigenvalue problem can be set up by expressing the radial variable y( 10 ) in the following 
time harmonic form: 

= y (r)e iut with wGl (3.50) 
where lo is the frequency of radial pulsations. With the introduction of the harmonic functional 



dependence d3.50t and the boundary conditions that later we discuss, the wave equation (13.441) 
becomes a Sturm-Liouville problem: 



l( P d f )-:-,,» + .=.r) f *-«. 



(3.51) 



where the squared frequency uj 2 is a free parameter. The solutions of this linear ODE form a 
countable set of discrete eigenvalues J 1 . 

In order to solve numerically this equation, one can transform it to a set of two first order 
ordinary differential equations by defining the new variable z^ 1 * ) = Py^ 19311591 . 



,(i,o) 



,(i,o) 



Z (1,0) 

p 



'W + Q)y 



(i,o) 



(3.52) 
(3.53) 



The boundary conditions associated with these equations are given by the regularity of perturba- 
tive equations at the origin and the vanishing of the Lagrangian perturbation of pressure on the 
surface. At the origin the condition (13.351) for ^ 1 >°) leads to the following behaviour: 



/ ( 1 -°)=y r 3 + O(r 5 ), 



,(i,o) 



y (i,o) 



+ 0(r 3 



(3.54) 
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that with equation (13.521) provides the following expression: 



yo 



3P 



(3.55) 



On the other hand, the vanishing of the Lagrangian perturbation of pressure on the stellar surface 
leads to the expression d3.42t . and then to the following condition: 



r=R 



0. 



(3.56) 



Various numerical methods can be adopted to integrate the system of equations (I3.52t - fl3.53l l. As 
we will see in section (I6.3.3t . we use the "relaxation method" ll9"Tll . 



3.4 Polar non-radial perturbations 

Linear non-radial perturbations of a non-rotating barotropic star can be described in the GSGM 
formalism by a set of metric and fluid gauge-invariant perturbation fields: 

{k^\k^\a^,^\J^} , (3.57) 

where all these tensors live on the submanifold M 2 and their definitions can be found in sec- 
tion !2.2l The metric perturbations are expressed by the two quantities k^g and /c^ ' 1 ), where the 
2-symmetric tensor k&B can be decomposed in the frame basis {u A , n 4 } as in equation ( I2.81I ). 
by giving the following expression: 



/ ( x (0,l) +fc (0,l)) e 2<£, _^(0,l) e $+A 

*S? = ! • (3.58) 

y _^(o,i) e *+A ( x (o,i) + fc (o,i)) e 2A 

In the previous equation we have used the definition ( 12.821 ) for the gauge invariant scalar pertur- 
bation x'' ' 1 ^ and the validity of the Einstein equation (I2.83t . which, for I > 2, sets a null value for 
the metric perturbation r/ ' 1 ). On the other hand, the fluid motion is described by the two gauge- 
invariant variables -yC 3 ' 1 ) and c^ 0,1 ) defined in (I2.73t and <2.72t . which describe respectively the 
velocity perturbation along the radial coordinates r and the latitude of the star. The perturba- 
tion of the mass-energy density is instead described by the variable u/ 0,1 \ see equation (12.74b . 
The equations of non-radial perturbations assume a simpler form when the enthalpy perturbation 
replaces the density perturbation ITTII951 80]. The enthalpy perturbations is so defined: 



ff (o,i) = _^_ w (o,i) < (3.59) 
P + P 

The six metric and fluid perturbations 

x (o,i) fe (o,i) ^(0,1) ; 7 (o,i) ( a (o,i) ) H (o,i) {3 60) 



CHAPTER 3. LINEAR PERTURBATIONS OF COMPACT STARS 



40 



are not all independent variables. Linear non-radial oscillations in the interior spacetime have a 
dynamics that can be described with two degrees of freedom, as shown in (25 . 56]. The external 
spacetime instead is a perturbed Schwarzschild spacetime where the single degree of freedom can 
be described by the Zerilli function 11151 . In the frequency domain approach, Chandrasekhar and 
Ferrari [25| were able to describe stellar oscillations in terms of pure metric perturbations. How- 
ever, the fifht order system of equations, which they derived in the "diagonal gauge", contained 
also one spurious solution. By using the Regge-Wheeler gauge, Ipser and Price 1551 succeeded 
to determine a fourth order system of equations, which contains only metric perturbations, and 
to clarify the origin of the additional solution of the diagonal gauge, which can be eliminated by 
a non-trivial gauge transformation. These researches then suggested that the fundamental infor- 
mation of the non-radial perturbations of relativistic stars are contained in the dynamics of the 
spacetime, as in the case for the black hole perturbations. As a result, other authors investigated 
the time evolution of non-radial oscillations with the same strategy, i.e. by setting up a system of 
PDE for these two metric degrees of freedom [96l l9"%ll . 

The GSGM polar perturbative equations that are reported in section l2~2l have been determined 
following the same method. The two independent perturbations are given by the metric variables 
1^.(0,1)^ ^(0,1) | th a t satisfy the two coupled wave equations (12.841) and (12.85b respectively, which 
are sketched here: 



The other four quantities can be then derived by using some of the other seven Einstein and 
conservation equations d2.86t -( l2~^2l . 

An analysis of the wave equation characteristics in the expressions (I3.61t and (13.621) allows 
us to interpret the physical properties of these two quantities. The perturbation x*- ' 1 ^ describes 
the gravitational degree of freedom that propagates according to the wave equation in all the 
spacetime. The interpretation of the variable k^ ' 1 ^ is less clear as it is a metric perturbations 
that propagates in the stellar interior at the speed of sound. In the exterior, equation (I3.62t losses 
its wave propagation character and another equation is required to evolve the variable A^ 0,1 ), see 
reference [96|. This problem can be avoided by adopting outside the star the Zerilli formulation. 

The polar non-radial oscillations can be also studied with a different set of independent vari- 
ables. In the work of Allen et al. lHH and Nagar et al. |81 ] for instance, the enthalpy perturba- 
tions H^ 1 ' ' has been evolved together with the two independent metric perturbations in a system 
of three partial differential equations. This system is well defined inside and outside the star al- 
though the enthalpy perturbation vanishes in the exterior. In reference lIHTTl . the authors use the 
GSGM formalism to set up a system of two hyperbolic and one elliptic equation for the three vari- 
ables {x^ 0,1 \ k^ ' 1 *} , ii^ ' 1 )}. The two hyperbolic equations are the gravitational wave equation 
(I2.84t for the variable x^ ' 1 ^ an d the sound speed equation for the enthalpy H^'°\ The elliptic 
equation is the Hamiltonian constraint ( 12.921 . which is used to update the metric variable A;( 0,1 ). 
The choice of solving the Hamiltonian constraint for one of the variables allowed them to have 



k 




(3.61) 



(3.62) 
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more control on the numerical errors and stable long time evolutions. 

In our work we have chosen to adopt this approach, as we think that higher control of errors 
and stable evolutions are crucial for having an accurate analysis of the non-linear coupling be- 
tween the radial and non-radial oscillations. In fact as will explain more in details in section l4~2l 
we argue that with this system of equations we can reduce the numerical errors also for the evo- 
lution of the second order perturbations. In addition, longer simulation times will allows us to 
investigate more accurately the interaction between the first order perturbations. 

The three equations for the interior spacetime are then given by the following expressions: 
Gravitational wave equation: 



S 



(0,1) 



+ e 2(*-A) 5 (o,i) +e 2(*-A) 



(5d> ir - A, r ) + - ( + * 2 r + ±£) 



+ - ( $ r (5 + 4$ r r) + 3A, r + 



2- G(/ + l) + 2)e 2A \ 



5 (o,i) 



0. 



(3.63) 



sound wave equation: 



-H 



(0,1) 



2 

- + 2$ 

r 



+ 



1 



[l + 3c 2 s ) (A r + d> ir )-c2fc^e 2A 



- A, P j c z s 
#(o,i) _ 



t2 



rS (o,D 
0, 



- -2$ 2 r + [(3*, r + A, r ) r + 1 - e 2A ] - 
Hamiltonian constraint: 

k<$>V S^+(--A jr )k^ + ^(A >r + $ r )H^ + \[(l-l(l + l))e 
\r J ' rCg r z 1 

+2A r r - 1] k^ l) - — \l(l + l)e 2A + 4 - 4A r r) 5 (0 ' 1} = . (3.65) 
2r L ' J 



1.(0,1) 

(3.64) 



2A 



The exterior spacetime is a perturbed Schwarzschild solution. The equations (13.631) and (13.651) . 
where all the fluid quantities present in the background coefficients and the enthalpy perturbation 



vanish, remain well defined. On the other hand, the sound wave (I3.65t is not defined there and 
the system of equations reduces to the two following equations: 



_ S (0,l) +e2( $-A )s (0,l) + e 2* 



6M 



AM 



S<2>*> 



M 



2M 



1 



2M 



,2A 



+ 



1(1 + 1) 



,2A) k (0,l) 



5(0,1) 
(3.66) 



- 2A (k^ - S^) 



+ 



2 2M 1(1 + 1) 



2r 



5(0.1) =0. 



(3.67) 
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From the previous two equations we can deduce the existence of a single degree of freedom for 
the exterior spacetime. In fact, the two metric variables have to satisfy the Hamiltonian con- 



straint (I3.67t . Zerilli showed in reference [115] that the propagation of the gravitational wave can 
be described by a single wave equation, afterward known as the Zerilli equation: 

- + e 2(*-A) z (o,i) + ^ e 2* z (o,D _ VzZ (o,D = o . (3 . 68) 

The Zerilli function, which represents the only degree of freedom of the external spacetime, is a 
gauge-invariant quantity related as follows to the GSGM metric variables S^ 0,1 ), k^ 0,1 >: 



Z (o,i) 



4r 2 e -2A 



1(1 + 1) [(l + 2){l-l)r + 6M\ 



(3.69) 

The function Vz is the Zerilli potential: 

v - f i - — \ ni<<ai - 2 - )2rS + 6 ^ ni - 2 ) 2Mr " 2 ± 36 ( n ' ~ 2 ) M2r + 72m3 

Z ~ ~ V ~~) r 3 [(m -2)r + 6M] 2 ' ( ' 

where ni = 1(1 + 1). 

The energy radiated at infinity in gravitational waves can be determined though the Zerilli 
function [29] with the following equation: 

dE m_ i (z + 2)! m 
dt - 647r ^— ' (1-2)1 1 ' 

which is valid for I > 2. 

The boundary conditions for the polar non-radial perturbations will be described in sec- 
tion !4.2.2l as a particular case of the polar non-linear Ae perturbations. 



3.5 Axial non-radial perturbations 

The linear axial, non-radial perturbations of a non-rotating star are described by the following 
metric perturbations: 

ds 2 (°4) = hf 1} Sa(dx A dx a + dx a dx A ) + (S a:b + S a:b ) dx a dx b , (3.72) 
and from the only fluid perturbations existent in the axial sector, i.e. the velocity perturbation: 

5u^ = (0,(3^S a ) , (3.73) 

where the scalar function Z^ ' 1 ) is a gauge invariant variable (see section 0, and in the previous 
two equations we have used the axial basis of the tensor harmonics defined in section |2~ 2. 21 

In the GSGM formalism the dynamical information of the spacetime is completely described 
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by the gauge-invariant master function ^I^ 0,1 ) defined in equation (12.95b . that in a static back- 
ground is given by the following expression: 



$(04) 



r k 



i.(o.i) 



"1A 



k 



(o,i) 

0,r 



(0,1) 



-(*+A) 



(3.74) 



while the perturbation Z^ ' 1 ) accounts for the amount of differential rotation present in the star 11091 
The system of perturbative equations can be determined by introducing the static background 



quantities given in the expressions (!3.13l) . J3.14t and (I3.15t into equations d2.93l > and < I2.94I >: 



_$(o,i) + ^(0 1) 



. ,_ , 6M 1(1 + 1) 
4vr (p - p) + — 



=2*^(0,1) 



,(0,1) 



+ 167T r 
0. 



„2<I>+A 



(3.75) 
(3.76) 



where we have introduced the tortoise coordinate dr* = e*~ A dr, in order to formally simplify 
equation (I3.75t . Furthermore, we have re-defined the axial velocity /?( 0,1 ) with the following 
perturbation: 

^(o- 1 ) = (p + p) 0^ . (3.77) 

The introduction of this new function is motivated mainly by two purposes: i) the equation at 
linear and non-linear order are simpler, and ii) for a poly tropic equation of state, (5^^ vanishes 
on the surface of a static star. This property will be very useful later in the numerical integration 
of non-linear axial oscillations, where the axial perturbations appear in the source terms. The nu- 
merical integration seems cleaner and more reliable with this new variable. From equation ( 13.761 ) 
emerges the stationary character of the linear axial velocity perturbation /iK ' 1 ) and then of Z^ ' 1 ), 
thus the amount of stellar differential rotation is completely determined by the spatial profile of 
the initial condition. Furthermore, the time independence of /3( 0,1 ) allows us to divide the solu- 



tion of the master equation (13.751) in two parts: i) the true dynamical degree of freedom of the 
spacetime, namely the gravitational wave, and ii) a stationary solution which is related to the 
dragging of the inertial frame due to the differential stellar rotation. Thus mathematically the 
general solution is given by 



^(0,1) = ^(0,1) , ^(0,1) 

^ ham. V ' 



(3.78) 



where the propagation of the gravitational wave ^ j.^ is carried out by the homogeneous equation 



associated with the PDE (13.75b. which is: 



^(O: 1 ) + $(o,i) 

At ' , r*r* 



+ 



4-7T (p ■ 



. 6M 
P) + — 



,(o,i) 



1(1 + 1) 



,2*^(0,1) =Q 



(3.79) 



On the other hand, the stationary metric profile \&p ' can be determined by a particular time- 
independent solution of the master equation. In this case, the system of equations is given by an 
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ordinary second order equation for ^p 0,1 ^ and a trivial evolution equation j3 <y0 ' 1 ^: 



p,rr ' 



3(0,1) 



47T (p - p) r H 2~ 



p. r 



_* , 6M 1(1 + 1) 
An (p - p) + — 



e 2A^(0,l) 



+ 167T r 
0, 



e- 2A ^ + (^npr + ^Pj ft ' 1 



,3 A 



(3.80) 
(3.81) 



The particular solution vt^, ' 1 ^ is related to the component kg of the gauge invariant tensor 



through the expression (12.96ft . The expression that connects this metric variable with the frame 
dragging function u>(r,9) is the following [109]: 



r 2 sin 2 6uj(r, 9) = 5g 



t(j> 



El,lm aim 
K <-></> ■ 



(3.82) 



1 1 n 



The gauge-invariant harmonic component k l ™ coincides with the metric component h l ™ in the 
stationary configuration, and more generally for time-dependent cases when the Regge-Wheeler 
gauge is adopted. Alternatively, the solution k l ™ and the related frame dragging u> lm can be 
determined directly by the following ordinary differential equation: 



. AM 1(1 + 1) 
8tt(p + p) + — 3 — 



,(0,1) 



16ne* ft ^ , (3.83) 



which has been derived by introducing the definition (13.741) into the expression ( I2.96I ). 

The axial gravitational radiation in the exterior spacetime produces small perturbations of the 
Schwarzschild solution. The propagation of these small ripples of spacetime is studied with the 
solution of the Regge-Wheeler equation [93], which is the equation one obtains by adapting the 



system of equations d3.75t and (13.761) to the exterior, 



0, 



(3.84) 



where the Regge-Wheeler potential is: 



V, 



(RW) 



r ) 



1(1 + 1) 6M 



(3.85) 



and r* is the usual Regge-Wheeler tortoise coordinate r* = r + 2M ln[r/(2M) — 1]. 

The boundary problem for the axial linear perturbations is complete when we specify the 
boundary conditions at the origin, at the surface and at infinity. The requirement of regularity of 
perturbation fields and equations at the origin gives 



3(0.1) 
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(3.86) 
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The junction conditions on the surface lead to the continuity of metric variable ^, of its time 
derivatives, and of the following expression [74]: 



In the case of a barotropic equation of state, the pressure and mass-energy density vanish on the 
surface, and the condition d3.87t induces the continuity of . At infinity, we impose for equa- 
tion d3.79t the Sommerfeld outgoing boundary condition in order to isolate the physical system 
under consideration. In addition, since at infinity the effects of the dragging of the inertial frame 
disappear, we set for equations (13.801 and ( 13.831 ) a vanishing value for ^p ' 1 ^ and ko respectively. 



For I > 2, The odd-parity master function ^ is related to the emitted power in GWs at infinity 
as |82l 



where we have explicitly restored the (I, m) multipolar indices and we have indicated with an 
overdot the derivative with respect to the Schwarzschild coordinate time. 




(3.87) 




(3.88) 



Chapter 4 



Non-linear Oscillations of Compact 
Stars 

Non-linear oscillations of relativistic stars have recently attracted great interest. The achieved 
improvements in numerical relativity now enable us to simulate non-linear evolutions of various 
physical systems such as core collapse l33ll34l . non-linear oscillations of differentially rotating 
stars [104], non-linear saturation of r-modes fj"5ll7Tll72"l ll051 non-linear radial pulsations [102], 
etc. 

In perturbation theory, the second perturbative order has been developed for studying Schwarzschild 
black holes [.^SJED- In particular, the authors treated the polar and axial quadrupole perturba- 
tions and applied this framework to the study of symmetric and asymmetric collisions of two black 
holes. The perturbative results have been compared with those obtained in numerical relativity, 
showing an unforeseen accuracy even for non-linear regimes. 

Perturbative techniques for studying non-linear oscillations have been used also for studying 
linear radial and non-radial perturbations of a slowly rotating star. In this case the rotating config- 
uration of the star is treated perturbatively [52], where the perturbative parameter is given by the 
dimensionless ratio Q,/Q,k, where Q, is the uniform angular velocity measured by an observer at 
infinity and Q,k the Keplerian angular velocity that determines the mass shedding limit. Several 
works have been carried out in the slow rotation approximation and with numerical codes devel- 
oped in numerical relativity for determining the effects of the rotation on quasi-normal modes 
(QNM). The splitting of the non-axisymmetric QNM frequencies and the presence of gravita- 
tional wave instabilities, i.e. the Chandrasekhar-Friedman-Schutz instabilities 1 241 l39l are two of 
the most relevant effects due to the rotation of the star (see review [ 103]). 

Another interesting property of the non-linear harmonics is that their frequencies may come 
out as composition frequencies, i.e. as linear combinations of the linear mode frequencies 1 11041 
11141 l65ll . This aspect could be interesting if two linear modes belonging to different families 
have nearly or exactly the same frequencies. In this case the associated non-linear harmonics 
could emerge at lower frequencies than the respective linear modes and with an amplified signal 
due to resonance effects. As long as they reach a high efficiency, these simultaneous effects 
could produce non-linear gravitational radiation in the sensitivity window of the new generation 
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of Earth based laser interferometers and mass resonant antennas. In this chapter, we derive for 
the first time the perturbative equations for studying the coupling between the radial and non- 
radial perturbations of a non-rotating compact star, where the stellar matter is described by a 
perfect and barotropic fluid. In section 14.11 we describe the method used for getting the non- 
linear perturbations, i.e. the 2-parameter relativistic perturbation theory on the GSGM gauge 
invariant formalism. The perturbations that describe the coupling between radial and polar and 
axial non-radial perturbation are introduced in sections l4~2l and l4"31 respectivelv. where we discuss 
also their boundary conditions. 

The results obtained in these sections have been presented in a first paper for the polar pertur- 
bations lIS^ll. while the axial sector is the subject of a second paper [ 90 1 . 

4.1 Coupling between radial and non-radial stellar oscillations 

The dynamics of the non-linear oscillations that describe the coupling between the radial and 
non-radial perturbations can be studied with a system of equations similar to the linear non-radial 
perturbations. In fact as explained in section 12.31 at second perturbative order the equations 
preserve the differential operator of the first order equations, and in addition they have quadradic 
source terms made up of the first order perturbations. At this perturbative order, the first order 
perturbations are already known by evolving the perturbative equations introduced in the previous 
chapter, and then behave as sources that drive the Ae perturbations. 

For the determination of the solutions of these inhomogeneous partial differential equations 
the spherical symmetry of the radial pulsations is very helpful. In fact as we see later, the variables 
and equations of interest can be obtained by expanding the gauge invariant formalism of GSGM 
with the 2-parameter perturbation theory. This approach allows us to split the time-dependent 
and spherically symmetric spacetime of GSGM in a static background and a radially pulsating 
spacetime. The consequences of this splitting will then be worked out also for the non-radial 
perturbations. 

The GSGM is valid for first order non-radial perturbations on a general time-dependent and 
spherically symmetric spacetime that we denote as follows: 



where M. is a manifold and Tq is the set of geometrical and physical tensor fields defined on this 
spacetime. In this formalism, the perturbed spacetime is a continuous one-parameter family of 
spacetimes diffeomorphic to the physical one, 



M G = (M,T G ) 



(4.1) 



M e = (M,T e ) , 



(4.2) 



where the strength of the non-radial perturbations is controlled by the perturbative parameter e, 
and T e are tensors that describe the geometrical and physical properties and are defined on the 
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physical spacetime. Any tensor field can then be perturbatively expanded as follows: 

T = T(°) + eT« + 0(e 2 ) , (4.3) 

where 7~(°) are tensor fields which describe the quantities of the time-dependent and spherically 
symmetric spacetime, and 7~W their non-radial perturbations. 

In our approach, the main point is the identification of the time-dependent and spherically 
symmetric spacetime Mq with a radially pulsating spacetime. In particular, since we are inter- 
ested in pulsations of small amplitude around an equilibrium configuration given by a non-rotating 
star, we can treat Mq perturbatively. As a result, the tensor fields associated with the Mq space- 
time can be written in the following way: 

T (0) =T + AT (1 - 0) + 0(A 2 ). (4.4) 

where T and all the quantities with the upper bar represent variables of the background that we 
consider in this work, i.e. a static star. Consequently T^ 1,0 ) denotes the one-parameter radial 
perturbations of a static star, where the strength of radial perturbations is controlled by the pertur- 
bative parameter A. 

In accordance with the new interpretation of the GSGM background, we can distinguish into 
the physical spacetime of non-radial perturbations M e the quantities that are dependent and in- 
dependent on the radial perturbative parameter A. The independent part will be denoted with the 
symbol T^ 0,1 ) and represents non-radial perturbations on a static stellar background. The second 
part, which depends on the perturbative parameter A, describes the corrections to the non-radial 
perturbations T^ ' 1 ) due to the radial pulsations. As a result, it represents the non-linear perturba- 
tions which describe the coupling between the linear radial and non-radial perturbations that will 
be denoted with T^ 1 ' 1 '. Any tensor can be then expanded as follows: 

T (D =T (o,i) +Ar (i,D +0(A 2 ) , (45) 

Hence, when we introduce the expressions (14.41) and (14.51) into equation (14.31) we find the following 
2-parameter expansion: 

T = T + A + eT^ 1 ) + Ae + 0(A 2 , e 2 ) . (4.6) 

It is worthwhile to remark that according to the multi parameter perturbation theory, the pertur- 
bative fields in the expression (14.61) are tensors defined on the background spacetime. 

This strategy is very useful for saving calculations and setting up a boundary initial-value 
problem for the description of these non-linear perturbations. In fact, we can determine the non- 
linear perturbative fields and the related systems of equations by introducing into the GSGM 
objects the expansion in the second parameter A, as shown in equation (14.41) and (I4.5t . The 
desired quantities will then be selected by virtue of their perturbative order Ae. Therefore, this 
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approach is a shortcut that prevents us performing a 2-parameter expansion of the Einstein and 
conservation equations directly from the static background quantities. 

Another important aspect of this method is that we can identify the gauge invariant variables 
for the metric and fluid nonlinear coupling perturbations. The method will be described in detail 
in chapter[5] we can give here only a qualitative analysis. Let's assume that of equation d4.5t 
is a GSGM gauge invariant quantity. Thus, we can immediately deduce that also the variable 
7~(04) j s gauge invariant. In fact, it is a first order non-radial perturbation of a static background, 
which can be treated as a subcase of the GSGM formalism. Let's now focus our attention on 
the Ae perturbation T^ 1 ' 1 '. If the gauge for first order radial perturbations is not fixed, T^ 1 ' 1 ) 
is not in general gauge invariant at order Ae. On the other hand, for a fixed radial perturbation 
gauge, equation (I4.5t may suggest that the gauge invariance of T^ 1 ' 1 ' can be derived by the gauge 
mvanance of TW and T^ 1 ). In chapter^ we will prove this property in detail by studying the 
structure of T^ 1,l \ which arises from the 2-parameter expansion of the GSGM formalism and by 
using the Ae gauge transformations for a fixed radial perturbation gauge. 

In the last part of this section, we will describe an alternative procedure for studying the 
non-radial perturbations on a radially pulsating star and will compare it with the 2-parameter 
perturbation theory. In this alternative approach, the GSGM time dependent background can be 
still treated perturbatively as in equation ( 14.41 ). whereas the quantity 7"M, i.e. a first order non- 
radial perturbation of a radially pulsating background, is now considered as an entire term without 
performing any further perturbative expansion as in equation d4.5t . The structure of the Einstein 
equations is then: 



^ NR. 



7~W 



+ \E 



r (i,o) g, T (i) 



+ 0(A 2 ) =0, 



(4.7) 



where we have not explicitly written the non-radial perturbative parameter e, and T represents 
both metric and fluid variables. We can notice that a part of equation (14.71) is governed by L atr, 
which is the linear differential operator d2.103t that acts on the first order non-radial perturba- 
tions (I2.106I) . while E is a linear differential operator that describes the remaining part of equa- 
tion ( I4.7I ). In addition, it is worth noticing that the error O (A 2 ) in equation d4.7t is given exclu- 
sively by the expansion of the time dependent background (I4.4I) . 

The perturbative equations ( 14.71 ) are now a system of homogeneous partial differential equa- 
tions, where T^ 1 ' is the unknown and the radial pulsating part of the time dependent background 
appears in the second term. Furthermore, it is worth noticing that the radial perturbative pa- 
rameter A explicitly appears in the equation and controls the strength of the radial pulsations. For 
A = 0, equation d4.7t describes first order non-radial perturbations on a static background d2.106t . 
In order to understand whether this equation is equivalent to the coupling equation d2.108t . we 
introduce the expansion < !4.5t in equation d4.7t obtaining: 



Lnr 



+ E 



+ X 2 E 



r (i,o) T (i,i) 



T (i,o) g, T (o,i) j j 
+ O (A 2 ) =0. 



(4.8) 
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At first order in e we have equation (12. 106b : 

l nr \tW] =0, 



(4.9) 



that describes the first order non-radial perturbation on a static star. At Ae order we get the 
perturbative equations d2.108t : 



Lnr 



+ E 



T (i,o) g, r (o,i) 



0. 



(4.10) 



which describe the coupling of first order radial and non-radial perturbations on a static star. In 
addition, a third part of higher perturbative order A 2 e is implicitly contained in equation (I4.7t . i.e. 



E 



T (i,o) g, T (i,i) 



0. 



(4.11) 



If we neglect terms of order eA 2 the solutions of equations (I4.7t and (I4.10t are equivalent. How- 
ever, we think that the 2-parameter perturbation theory can give cleaner results. First, it directly 
gets rid of higher order information, as for instance the extra terms of order A 2 in equation d4.8l >. 
Even though this term is of order A 2 e, it is formally contained in equation (14.71 . Second, it enables 
us to work with inhomogeneus partial differential equations (12. 1086 . where the linear differential 
operator is defined on the static background. On the other hand, the differential operator of the 
homogeneous equation d4.7t has a part which is governed by radially oscillating quantities (see 
second term of equation \4.1V ). Therefore, we expect that the simpler differential operator of 
equation ( 14.101 ) will produce less numerical problems during the implementation of the code. 

4.1.1 GSGM formalism on a radially oscillating star 

In order to implement the GSGM formalism one has to specify all the quantities A2.34I ). ( 12.351 ). 
( 12.391 ) and (I2.41t that describe the time-dependent and spherically symmetric spacetime Ma- 
in our approach this spacetime is considered as a radially pulsating spacetime which is treated 
perturbatively. Therefore, all the variables defined on it must be expanded in a static and radial 
pulsating part as in the equation d4.4t . 

The frame vector field basis \u^' a , A } of the submanifold Mq assumes then the follow- 
ing expressions: 



ii 



(0)A 



1_ A L(i.o)_* 



(1,0) 



n 



(0)A 



Ae-V' 0) ,fl-A^ ]e- A 



(4.12) 
(4.13) 



These two vectors satisfy the ortho-normalization up to A 2 perturbative order: 



(0)A (0) 



U v ' u 
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(0)A (0) 



n v ' n 
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n v X'=0(X z ). (4.14) 
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With these definitions the metric and the completely antisymmetric tensors of the spacetime Mq 
take the following form: 



(0) 

9ab 



V°U°> +n (0) n (0) 
U A U B + n A n B » 



,(°) 
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(4.15) 



The action of the frame derivatives on a scalar perturbation /W = + Xf^ 1 ' 1 ' defined on the 
radially oscillating star Mq is defined by the following expressions: 
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X 



(1,0) 



- A /, ( r ' 1} +0(A 2 )- 



(4.16) 



(4.17) 



Finally, we consider the remaining scalar functions defined in this formalism, which are the com- 
ponents of the vector va = r'/r with respect to the frame basis, i.e. U, W, and the divergence 
of the two basis vectors, i.e. fj,, v. The expansion of the manifold Mq leads to the following 
expressions: 



U<P) 
,(°) 



t (0)^(0) =A £_ 7 (l j0 ) +o(A2)) 

n (o)A v (o) = (l_ x ^_) e -A + (A 2 ), 



/./> 



MA 



2r 



A ( 7 (1 ' 0) e" A ) ^ + 0(A 2 ) , 



,(o) 



nf/ = $ r e~ A + A \ e-^l fi) + e" A 



„(1.0) _ ly(l.O) 
V,r 2 



1 



#, rX a.o) 



(4.18) 

(4.19) 

(4.20) 

+ 0(A 2 ), 
(4.21) 



where in order to simplify the function (I4.20t we have used the radial perturbative equation (I3.26t . 



4.2 Coupling of radial and polar non-radial perturbations 

The non-linear perturbations and the perturbative equations at order Ae are determined with the 
2-parameter expansion of the GSGM perturbation fields. The explicit expressions of all the metric 
and fluid perturbations will be given in chapter |5j where we address the gauge invariance issues. 
In this section, we focus only on the dynamical variables formed by the two metric scalar fields 
X^ 1 ' 1 ) and fc^ 1 ' 1 ) and the enthalpy H^- 1,1 ' . In this way, we can set up a system of perturbative 
equations with the same differential operator as the first order equations (I3.63t - (l3.65t . 
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In the GSGM formalism the even-parity metric perturbations are described by the 2-rank sym- 
metric tensor kP B and the scalar A^ 1 ) that are respectively defined in equations d2.59l > and (12.601) 



on the radially oscillating spacetime Mq- The symmetric tensor k AB can be then decomposed 
with respect to the vector frame |^^' n A^} °f ^G?> an d three gauge-invariant scalar functions 
rp-\cjp-) and if)^' can be defined (I2.81I) . The explicit expression of the gauge-invariant tensor 



^AB components can be determined by introducing the expansion of the basis vectors ( 14.121 ) 



and (147131) 



1.(1.1) 
v oo 

u(Ll) 

m 



x (i.D + + (2^0) _ x d,o)) ^(o,D + fc (o,i)) + 2 7 (i,o)^(o,i) 



^(i,i) + 2 7 d.o) r x (o,i) + k m\ + ^(1.0)^,(0,1) 
x (i,i) + fc (i,i) + x (i.o) Uo,i) + k (o,i)\ + 2 7 (i.o)^(o,i) 



e 2 *, (4.22) 



(4.23) 
(4.24) 



where we have used the definition ( 12.821) for the x perturbation 



(4.25) 



and the Einstein's equation ( I2.83t : (7/ 1 ) = 0). 

The other gauge invariant quantity fcW assumes instead the following form at Ae order: 



^(1,1) = K (1,D 



-2A 



p (i,D_ x (i,o) p (o,i) 



(4-26) 



where we have performed the further perturbative expansion in A of the vector 



(4.27) 



which according to its definition (12.651) is given by: 



(0,1) 

Pa 

(i.i) 
Pa 



,(o,i) 



h 



(i,i) 



r r (0,l) 

r<1 r (l,l) 
~~2 \ A 



(4.28) 
(4.29) 



4.2.1 Perturbative equations for polar perturbations 

The explicit form of the perturbative equations that describe the dynamical properties of the non- 



linear coupling can be derived from equations J2.84t -( I2T92| > by introducing the quantities (14. 121) 



(14.211) of the radially oscillating spacetime Mg and the perturbative fields, which are further 
decomposed as in equation (14.5b - The desired equations are given by the perturbative part as- 
sociated with the perturbative parameter Ae, and accordingly the equation (12.1 16t will have the 
following structure for any harmonic index (l,m): 



j NR 



B (M) w/i,i) 



s 



r(l,0) 
'00 



(0,1) 
Im 



Ml > 1 



(4.30) 
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This particular structure of the (1,1) equations is quite convenient in order to build a boundary 
initial-value problem and solve it numerically by using time-domain methods. The basic idea is 
that given a numerical algorithm capable of evolving linear non-radial perturbations, we can build 
an algorithm for our (1,1) perturbations by just adding source terms to the original algorithm. 
The time evolution of non-radial perturbations of a static star has been successfully analyzed by 
numerically integrating different systems of perturbation equations fTTl [%1 l8T1l . However, for 
the main features of our formulation, the scheme introduced by Nagar et al. [ 8 1 , 80 ] seems to be 
more appropriate for the implementation of a numerical code. In their scheme the Hamiltonian 
constraint is not just an error estimator for the evolution equations, as is usually done in many free 
evolution schemes. Instead, it is part of the system of equations and is solved at every time step 



for the perturbative quantity k, equation (12.60b . This provides some control of the errors induced 
by constraint violation. As a consequence, the resulting numerical code [81] is able to evolve 
non-radial perturbations for long times and is capable of estimating the damping time and mode 
frequencies with an accuracy comparable to frequency domain calculations. 

Therefore, since the system of perturbative equations for the linear and non-linear, non-radial 
perturbations have the same differential structure we can expect to control the numerical errors 
by solving the Hamiltonian constraint both at first and second perturbative order. Otherwise, if 
we do not use this hyperbolic-elliptic system of equations the errors accumulated from constraint 
violation would be double that in a standard computation of non-radial perturbations. Therefore, 
we expect that this scheme allows us to obtain accurate long term evolutions. 

In the stellar interior, we evolve a hyperbolic-elliptic system of three partial differential equa- 
tions for the three second order variables: 

■ (4-31) 

The two metric variables x^ 1 ' 1 ^ an d fc^ 1 ' 1 ) defined above, obey respectively a wave-like equation 
that describes the propagation of the gravitational radiation and the Hamiltonian constraint. The 
Hamiltonian constraint is an elliptic equation and is solved at any time-step to update the value 
of k^ 1 ' 1 '. The third variable is the second order fluid perturbations H^ 1 ' 1 ', whose expression can 
be derived by the perturbative expansion of the analogous quantity on the radially oscillating 
spacetime 

« m -M0> J1) - ,432) 

This procedure leads to the following expression: 



#(1,1) = ^^(1,1) 
P+P 



dp P + P 



0,(1,0)^(0,1) f (4 _3 3) 

P + P 



where we have introduced the perturbative expansion of the sound speed on the radially pulsating 
spacetime: 

c 2(0) =d 2 X d^ 6p (l,0) (434) 

dp 
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In particular gauges, the Regge-Wheeler one for instance, the gauge-invariant quantity a/ 1 ) co- 
incides with the gauge dependent perturbation cl^ 1 ) [see equation (12.741) 1 . and if W describes the 
enthalpy perturbation, 

5pW 



(4.35) 



p(o) + p (o) ' 

where 5p^ is given by the definition < I2.7 II . 

The other metric (ip^'V) and fluid (j^- 1 ' 1 ', a^ 1 ' 1 ') perturbations can be successively derived 
by solving the perturbative equations (12.861) . ( I2.91t and ( 12.871 . 

The wave equation for the variables x^ 1 an d tne Hamiltonian constraint for k^ 1 ' 1 ' are given 



by equations (12.841 and (12.921) respectively. On the other hand, the sound wave equation for 
the perturbations H^' 1 ' must be determined by using the perturbation of the conservation equa- 
tions given in l74"ll . To this end, we prefer first to operate within the GSGM framework and 
find the equation for the perturbation ffW on the radially pulsating spacetime Mq, and sec- 
ondly carry out the 2-parameter expansion and derive the equation for the non-linear variable 
H^K This equation is obtained as a linear combination of the time frame derivative of equa- 
tion ( 12.891 and the spatial frame derivative of equation ( 12.901 . After having introduced equations 
(12.861 . d2. 881 . (12.901 . (12.851 and (12.921 to reduce the number of perturbative unknowns and the 
transformation (14.321 . we have the following wave equation 



(4.36) 



where Th contains all the remaining terms (with derivatives of lower order). The complete equa- 
tion has been written in Appendix |E] It is worthwhile to remark that the wave equation ( 14.361 is 
valid in the GSGM framework for barotropic non-radial perturbations on a time dependent back- 
ground. In case of a static background, given the introduction of the static quantities (I3.13t - d3.15l l. 
it reduces to the equation used in the literature (see references llTTll96l [81 1). 

We can now write the perturbative equations for the stellar interior. We consider instead of the 
perturbative quantity x > which diverges like r as we approach spatial infinity, the perturbation 



variable 



X 



(1: 



x ) jr which of course is well behaved at infinity. This quantity satisfies the 



following gravitational wave equation: 



S 



(1,1) 



+ 



+ 



2($-A) 0(1,1) , JJ(*-A) 

1 



(5$ r - A, r ) Si 1 ' 1 ) + 



,2A 



+ $,r + 



r 



* r (5 + 4* „) + 3A r + *-('(' + !) + S (U. 

r / 



„2$ 



r 
S S , 



k^ 
(4.37) 



where Ss denotes the source term for this wave equation. In particular, the source term in the 
gravitational-wave equation, Ss, has the following form 



S S 



2e A-*^i) N ) + a 6 /fe(°' 1 ) 
(4.38) 
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where the coefficients ctj are just linear combinations of radial perturbations with coefficients 
constructed from background quantities. Their explicit form is given in Appendix ICl 

The perturbative fluid variable H^' 1 ' also satisfies a wave equation, but with a different prop- 
agation speed. We call this equation the sound wave equation. It has the following form: 



- + 2$ - - A , c« - $ 



1 

+ - 

r 



[1 + 3c 2 ) (A, r + *, P ) 



,1(1 + 1) 



+ 



-2$ 2 r + [(3$ r + A, r ) r + 1 - e 2 
and the source term can written as 



S 2A 

-s 



1 ^ 2 



(rS(i,i)) _ fc (i,D 



r5< 1 ' 1 )+^ 1 ' 1 ))j = e 2 *5 Hj 



(4.39) 



5 H = ^ijCO-i) + ^ (o i) + hH i°,D + 64jff (o,i) + &5iJ (o,i) + 6gA; (o,i) + ^ 



r(0,l) 



(0,1) 



7 (0,1) 



+ h 



k m _ ( r5 (o,D) 1 + bg ( r5 (o,D + fe (o,i)) + 5lo7 (o,i) + fen7 (o,i) 



+ W (0.1) +6l 3^(0,l) +6l4Q (0,l) j 



(4.40) 



where the coefficients 6j have the same structure as the a-i coefficients in (14.38b . Their explicit 
expressions can be found in Appendix ICl 

For the last perturbative variable, the metric perturbation fe( 1 ' 1 ) , we will use the Hamiltonian 
constraint instead of an evolution equation. After some calculations we get: 

fc(l,D _ 5 (1,D + (1 _ A ) fc (l,l) + _L (A + $ } ^(1,1) _ 1 r /(/ + l)e 2A + 4 _ 4A 1 

+ 1(1 + 1)) e 2A + 2K r r - l] = S 



(4.41) 



where Sjj am n is the source term for the Hamiltonian constraint. As in the previous equations the 
precise form of S H amii is: 

(fcJW - SfV) + c 2 k^ + C3k fV + C4 5(o,D + C5fc (o,D + C6i j(o,D + ^(0,1) 
+ csV' 10 '^ + c 9 7 (0 ' 1) . (4.42) 

The coefficients c%, in the same way as the coefficients a-i and 6, only contain radial perturbations 
^(i.o) anc j quantities associated with the static background. They are also given in Appendix lO 
It is worth remarking that the polar non-radial perturbation equations on a static background are 



obtained from equations (I4.37I) . ( 14.391 ) and d4.41t by discarding the source terms and replacing 
all the (1,1) perturbations with the corresponding non-radial (0, 1) terms. The sources are de- 



termined from first order perturbations. The radial perturbations from equations (I3.24I) - (I3.27I) . 
and the non-radial perturbations (described by the quantities 5( 0,1 ) , A^ 0,1 ) , and H^ ' 1 ^) from the 
first order analogous of the above system (see reference IHT1 ). and the equations A2.86I ), ( 12.911 ) 
and (12.871) adapted to a static background to get i/^ 0,1 ), 7 (°> 1 ) and c/ 0,1 ). 
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The stellar exterior is described by a Schwarzschild spacetime on which gravitational waves 
carry away some energy of the stellar oscillations. All fluid perturbations vanish outside the star 
and the radial perturbations do the same because of Birkhoff 's theorem. Therefore, the source 
terms in our perturbation equations disappear. Only the metric perturbations survive, and they 
satisfy the gravitational wave equation d4.37t and the Hamiltonian constraint (14.41b . which take 
the following form: 



c(M), 2(<£>-A)~(l,l) , p 2<K 
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6M^(i,D 
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(4.43) 



-2A 



_ 5(1.1) 



+ 



1(1 + 1) 



+ + 5(1,1) =0 . 
2r 



(4.44) 



It is worth mentioning that the above equations coincide with the equations for non-radial pertur- 
bations of a static stellar background outside the star, as expected. 

On the other hand, Zerilli showed that the even-parity perturbations of a Schwarzschild back- 
ground have just one degree of freedom, and therefore can be described by just one variable, the 
Zerilli function, satisfying a wave equation. At order (1,1) the Zerilli function can be built from 
the two metric perturbations S^ 1 ' 1 ) and fe^ 1 ' 1 ) and their derivatives, as at first order an d 
is given by 



4r z e 



2„-2A 



1(1 + 1) [(l + 2)(l-l)r + 6M] 
It satisfies the Zerilli equation II 16111 15ll 

,(1,1) 
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z (i,D + e 2(*-A) Z (i,D + ™ e ™ Z (y) - v z z^ = , 



(4.45) 



(4.46) 



where Vz is the Zerilli potentiaK I3.70t . 

Finally, we can determine the power of the gravitational radiation emission at infinity by using 
the following expression [ 29 ] 



(It ~ 2^ 



a + 2)!, ^(i,i),2 



64tt ^ (Z - 2) ! 1 Zm 

l ,m 



(4.47) 



for Z > 2. 



4.2.2 Boundary conditions for polar perturbations 

In this Section we discuss the boundary conditions at the origin, infinity and at the stellar surface 
for the Ae perturbations describing the coupling of radial and polar non-radial modes. With re- 
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gard to the outer boundary, we locate it far enough from the star and we impose the well-known 
Sommerfeld outgoing boundary conditions on our perturbative fields. 

At the origin, the boundary conditions are just regularity conditions on the perturbative fields, 
which can be obtained by a careful analysis of the equations that they satisfy. The analysis of 
Taylor expansions of the differential operators that appear in our equations near the origin leads 
to the following behaviour for the non-radial perturbations [48 1: 

I > 0, S w ~ r l+1 fcW ~ r l ^W~r m , (4.48) 
I > 1 , 7 (1) ~ r 1 ' 1 ~ r l a ~ r l , (4.49) 

where the upper index (1) in the previous quantities denotes the linear e and non-linear Ae per- 
turbations as in equation d4.5b . The conditions (I4.48b - (l4.49ft and the expressions d3.32l i- (l3.35t for 
the radial perturbations leads also to the regularity of the source terms. 

The stellar surface in spherically symmetric spacetime is a 1 -dimensional manifold embedded 
in the 2-dimensional manifold M 2 C M ( section l2.2.1t . In order to prevent S discontinuities on 
the energy momentum tensor the two fundamental forms, i.e. the induced metric tensor and the 
extrinsic curvature, have to be continuous on this hypersurface. In particular, when the perturba- 
tive approach is used, these continuity conditions have to be imposed at any perturbative order 
considered. A boundary condition must be imposed at the surface also for the matter variable H , 
which vanishes outside the star. 

Let £ be the surface of the static unperturbed star (i.e. r = R s ). The surface of the perturbed star 
can then be defined as 

S = f x (t) = x + A£ (1,0) + e^ ' 1 ) + Ae£ {1,1) : x £ e| , (4.50) 

where is a vector field that denotes the Lagrangian displacement of a fluid element due to 
the action of perturbations of order A physical requirement that follows from matching 

conditions is the vanishing of the unperturbed pressure p at the unperturbed surface E. In the 
same way, the corresponding boundary condition for the perturbed spacetime is the vanishing 
of the total pressure p = p + Xdp^ 1 ' ^ + edp^^ + Xedp^ 1 ' 1 ^ at the perturbed surface E. This 
condition turns out to be equivalent to the vanishing of the Lagrangian pressure perturbations on 
E, the unperturbed surface, at every order. The Lagrangian pressure perturbations are given by: 

6pW + £ km p, (4.51) 
Sp^ + £ koi} p, (4.52) 

+ (^(M) + \ {%,0) . ^,0.!)}) P+£^M m + ^oM ^ 

+ (^ (M ) " \ {%,o) • ^,o.d}) P. ( 4 -53) 



Ap (i,o) = 
Ap (o,i) = 
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where 5 and A denote the Eulerian and Lagrangian perturbations respectively, and we have used 
the lower order boundary conditions Apt 1 ' } = Apt ' 1 } = in order to simplify the condition 



We can then conclude that the boundary conditions for the fluid perturbations are described 
by the set of expressions given in d4.51t -( R3^ l. However, in practice, in many applications of 
first order perturbation theory, dynamical boundary conditions either for density or enthalpy per- 
turbations have been considered. This alternative boundary condition follow from the analysis of 
the time derivative of the condition (I4.52t (see [ 1 1 ] for more details). In our current development 
of the numerical implementation of the perturbative equations we are considering both types of 
boundary conditions with the perspective of analyzing which type works best for our formulation. 

Finally, the junction conditions for the metric perturbations can be determined by imposing 
continuity of first and second fundamental differential forms and their perturbations at the sur- 
face imi69ll86ll57ll. The explicit form of these conditions has been presented in [42, 74] for first 
order perturbations of a time-dependent stellar background. According to our interpretation of the 
GSGM formalism (see section l4~lT) the time dependent and spherical spacetime is identified with 
a radially pulsating spacetime, which is treated perturbatively. Therefore, in order to determine 
the correct expressions of the junction conditions for linear e and non-linear Ae perturbations on 
the static surface we have to perform an expansion of the GSGM matching relations [74] similar 
to that describe in equations (I4.52t -( l4.53t for the pressure. However, since we are carrying out 
the analysis in an Eulerian gauge, the implementation of the junction conditions for the non-linear 
perturbations requires some approximations 1 1 H 1 1 02ll . In fact, due to the movement of the stel- 
lar surface some perturbations can take unphysical values near the surface during the contraction 
phases of the star. This is the case for instance of the total density p, which can become negative. 
Furthermore, the low densities which are present at the outermost layers of the star can produce 
also some numerical errors in the simulations IBTll . These problems can be avoided by imposing 
the matching conditions for the non-linear perturbations not on the static surface, as we do for the 
linear perturbations, but slightly inside [ 10l||^. This approximation corresponds to neglecting 
less than one percent of the stellar mass, which does not produce significant changes in the wave 
forms and spectra of the gravitational signal. 

For first order polar non-radial perturbations on a static star, the junction conditions can be 
determined at the static surface r = R s as a subcase of the expressions given in reference f74ll . 
These relations provide the continuity of the following I > 2 scalar fields: 

fc (o,D 5 5 (o,i) > ^(o,i) 5 fc(o,i)' ; 5 (o,i)' ; (4 . 54) 

where we have also used the vanishing of the static pressure, density and speed of sound at the 
stellar surface. 

At second perturbative order, the non-linear perturbations are matched on the following hy- 
persurface that during the evolution is always inside the star, 



Sj c ={r = R jc : R jc < x(t)} , 



(4.55) 
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where x(t) describes the position of the perturbed surface as defined in equation (14.50b . Therefore, 
we can perform a 2-parameter expansion of the GSGM junction conditions [74] at S JC as well as 
we have done for the pressure in equations d4.52l i- J4.53l l. By using the continuity conditions for 
the first order perturbations on the surface S JC , we obtain the following continuity conditions at 
Ae non-linear order: 



5(1,1) 



^ , e(- A ¥^ + ^e A H^ , SW) . (4.56) 



Alternatively, one may use the "extraction formulas" [74] that relate the Zerilli function with 
metric perturbations at the stellar boundary. For linear polar perturbations on a static star this 
relation reads for I > 2: 



£(0,1) = rk (p,i) 



+ 



2r 4 



(1 + 2) (Z-l)r + 6m 



-2A 



.5(0,1) + fc (0,l)) 



-2A 



! 1.(0,1) 



. (4.57) 



The extraction formula for the non-linear Ae perturbations can be imposed on the hypersurface 
S JC . By using the first order continuity conditions we obtain: 



Z 



(1,1) 



+ 



2r 4 



(I + 2) (I - 1) r + 6m 



„ -2A , . e -2A 

^5(1,1) + k d,D) _ e _( k {y) 



+ ^ e A^(l,D 



(4.58) 



4.3 Coupling of radial and axial non-radial perturbations 

In this section we derive the perturbative equations that govern the non-linear coupling between 
first order radial and axial non-radial oscillations. Following the same method of the polar case 
in section 14.21 we perform a 2-parameter expansion of the perturbation fields in the odd-parity 



perturbative equations ( 12.931 ) and ( 12.94b . 

The odd-parity metric perturbation of order (1,1) can be expanded in tensor spherical har- 
monics as 

ds 2 ^ = h i l' 1) S a {dx A dx a + dx a dx A ) + (S a:b + S a:b ) dx a dx b . (4.59) 

while the axial velocity perturbation of the fluid is given by 

5u^ = {o,p^S a ) . (4.60) 

The information of the dynamical properties of the spacetime is completely contained even at this 
perturbative order in a scalar function \E'( 1 > 1 ). This axial master function can be determined by 
expanding equation (I2.95t as we showed for the metric in equation (14.5b . leading to the following 
expression at Ae order: 



*^ = [r(jfc^-^) 



+ 2k 



(1,1) 



-(»+A) _ (1,0)^(0,1) _ 



(4.61) 



CHAPTER 4. NON-LINEAR OSCILLATIONS OF COMPACT STARS 



60 



In the previous equation the Ae perturbative components of the gauge invariant vector fcy appear, 
see equation ( I5.19t . 

The system of perturbative equations is formed by the master wave equation ( 12.931 ) and by the 
conservation equation j2.94l >. which assumes the following form after a 2-parameter expansion: 



(i,i) 



+ * 



+ 167T r 



+ 



■ 6m 1(1 + 1) 



M 



e 2*^(l,l) 



e~ 2A ^ + Uirpr + ^)pW 



2<I>+A 



V 



(4.62) 
(4.63) 



where the explicit expressions of the source terms £^ and S/j are given in Appendix ICl 

In the previous perturbative equations, the axial velocity perturbation fyl^ has been replaced, 
as in the case of the first order perturbations, with the variable fit 1 ' 1 '. In this way we have a uni- 
form set of perturbative variables for any perturbative order. In order to get its explicit expression 
we can first define this variable for a generic time-dependent spacetime: 



/3 ( 



i) 



(4.64) 



and then introduce the perturbative expansion in two parameters. At order e the definition (I3.77t 
valid for linear axial velocity is restored, while at non-linear order, i.e. Ae, the following expres- 
sion is determined: 



l + c 



#(1.0) 0(0,1 



(4.65) 



It is worth noticing that due to Birkhoff 's Theorem, the source terms are present only in the 
stellar interior. Hence, the exterior is a Schwarzschild space-time perturbed by the gravitational 
waves, which can still be described by the Regge-Wheeler equation of Ae order: 



Consequently, the emitted power in GWs at infinity reads: 

1(1 + 1) 



(4.66) 



dE^ _ _L_ 
dt ~ 16vr 



E 

Lm 



(l-l) (1 + 2) 



(i,i) 

lm 



(4.67) 



where we have explicitly written the multipolar indices (I, m) with I > 2 and we have indicated 
with an overdot the derivative with respect to Schwarzschild coordinate time. 
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4.3.1 Boundary conditions for axial perturbations 

We still need to impose the correct boundary conditions to Eqs. ( I4.62t and (I4.63t at the origin, at 
the stellar surface and at infinity. At the origin a regularity analysis shows that 

£(1.1) ~ r l+1 ^r(i>i) ~ r l+1 . (4.68) 

The junction conditions at the stellar surface for the non-linear axial perturbations can be de- 
termined by using the GSGM equations for a time dependent and spherically symmetric space- 
time [44, 74], as we have already illustrated for the polar case in section 14.2.21 In the axial case 
the movement of the surface is due uniquely to the radial perturbations, as the axial perturbations 
at e and Ae perturbative orders can only produce a fluid motion in the tangential direction. The 
continuity of the induced metric, extrinsic curvature and their perturbations leads to the follow- 
ing results: i) continuity of xfj^ and its time derivatives (ip^)' , and ii) continuity at the stellar 
surface of the following expression: 

(V 3 #(!)) + 16vrr~ 2 /3 (1) . (4.69) 

According to our interpretation of the GSGM formalism (see section |4~TI . the correct expressions 
for the junction conditions on a hypersurface S JC placed inside the star can be determined with the 
procedure used in section l4~2.2l At first perturbative order in e, we get the same conditions for the 
linear axial perturbations expressed in section 1X51 For the non-linear coupling, the continuity of 
the master non-linear function vp(i>i) and its time derivatives must be imposed at S JC . In addition, 
the condition ( 14. 691 ) leads to the continuity of the following quantity 

e -A tf J*' 15 + 16vr r /3 (1,1) , (4.70) 

where we have used the continuity of \E'( 1 > 1 ) and the linear junction conditions. At infinity, we 
impose the outgoing Sommerfeld boundary condition. 



Chapter 5 



Gauge Invariance of Non-linear 
Perturbations 

This chapter is dedicated to the gauge invariance of non-linear perturbations of relativistic stars. 
The discussion is specialized to a particular class of second order perturbations, namely those 
that describe the coupling between the radial and non-radial oscillations. The construction of 
gauge-invariant quantities for these non-linear perturbations will be based on the same strategy 
that we have used for the determination of the perturbative equations, i.e. a 2-parameter pertur- 
bative expansion of the GSGM formalism, see section |4~TI This approach will be very helpful for 
identifying and building the non-linear perturbations with gauge invariant character directly from 
the gauge-invariant quantities of the GSGM formalism. 

The gauge invariance of relativistic linear and non-linear perturbations has been investigated 
in many works ll06H20l[T9l[T0 0, 83 ], where transformation rules have been presented for studying 
one or multi-parameter perturbative problems. In addition, within the multi-parameter perturba- 
tive framework, a technique for the construction of gauge invariant, non-linear perturbative fields 
has been introduced by Nakamura [ 83 , 84 ]. In the literature, there are different ways to define the 
gauge invariance of non-linear perturbations [20|[3U|. A perturbative tensor field of n order can 
be defined as gauge invariant "at n perturbative order" or "up to n perturbative order". For the 
former definition a perturbative tensor field is invariant only for the gauge transformations relative 
to the n perturbative order considered. Therefore, all the gauges of the previous 1, .., n— 1 orders 
must be fixed. The second definition is more restrictive and requires that the perturbative field is 
gauge invariant at any order up to the desired perturbative order. Thus in this case, all the gauges 
from the first to the n th orders are completely arbitrary. 

The non-linear perturbations that we adopt for the analysis of the coupling between the radial 
and non-radial oscillations will be gauge invariant for gauge transformations associated with the e 
and Ae non-radial perturbations with the restriction that the first order gauge for the radial pertur- 
bations must be fixed. This restriction is due to the absence of a gauge invariant formulation for 
the radial perturbations. In fact, as we reported in section 13.31 the GSGM formalism also fails to 
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provide a gauge invariant description for this class of oscillations. However, their physical prop- 
erties have been well described in both Eulerian and Lagrangian gauges (59JE0I1. In particular, 
the latter is more appropriate for describing the regions near the stellar surface II 01 H . 

In section 15.11 we present the method to build the gauge invariant perturbations at Ae order. 
In section 15.1.11 we illustrate this method for the axial perturbations while the polar sector is 
discussed in section l5".1.2l This method has been presented in the paper [89] for the polar pertur- 
bations, while the axial case will be treated in [90]. 

5.1 Construction of gauge-invariant non-linear perturbations 

Gauge transformations and gauge invariance in 2-parameter perturbation theory have been stud- 
ied in reference fT9l fl001. In section IXTl we have reported the procedure based on the Baker- 
Campbell-Hausdorff (BCH) formula given in the work of Sopuerta et al. COO], where the gauge 
transformation rules for linear and non-linear perturbative fields have been derived in a more 
general way. Let T be a generic tensorial quantity and T^' its perturbation at AV order. 

The gauge transformations of linear radial and non-radial perturbations T0-> Q ) and T^ 0,1 > are 
respectively given by the formulae (12.231) and (I2.24I) . 

f(i,o) = T ^+£^ 0) T b , (5.1) 
r(o,D = T^ + £^ 1} T b , (5.2) 

where T b is the background value of the tensor field T, and the tilde denotes the perturbative fields 
transformed by the gauge transformation. The vector fields £(i 5 o) an d £(0,1) WQ respectively the 
generators of the gauge transformations for the first order radial and non-radial perturbations. 
At second Ae order, the tensor fields T transforms according to equation (12.271 : 

fW> = rfcD + £ k01) + £ klfi) T<M + (j? g(M) + {£ kl0) , £ ({01) }) % , (5.3) 

where { , } stands for the anti-commutator {a, b} = a b + b a, and £(1,1) is the generator of the Ae 
gauge transformations. In the transformation rule (15.31) . the presence of non-linear terms, which 
contain perturbative fields and gauge transformation generators of the previous order, make really 
hard the construction of non-linear perturbations that are completely "gauge invariant up to the 
second order". A method to reach this purpose has been presented by Nakamura in reference 1531 . 
but here, we do not follow his approach as the calculations will be too laborious for the aim 
of this work. Furthermore, as long as a gauge invariant description for the radial pulsations is 
unknown we have to study these linear perturbations by choosing a gauge, for instance as we did 
in section 13.31 Therefore, we can derive non-linear perturbative fields which are invariant under 
a sub-class of second order gauge transformations, namely those where the first order gauge for 
radial perturbations is fixed. With this assumption the general gauge transformation d5.3t reduces 
to this simpler expression: 



f(U) = T^) + ^ ([U) T^) + £, (M) T fe , 



(5.4) 
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which is determined by setting to zero the radial generator £(i,o) in the expression (I5.3t . 

Now, we are going to explore the structure of the non-linear perturbations that we have used in 
chapter|4]and as these perturbations change under the gauge transformation d5.4b . The main idea 
behind our method is to build the Ae gauge-invariant variables starting from the gauge-invariant 
quantities of the GSGM formalism. 

Let C/M be any of the gauge-invariant metric or fluid quantities in equations d2.59t - d2.64l) 
and d2.72t -( l2~74l for the polar sector or in equations d2.66t -( l2~6^l and ( 12.701 ) for the axial case. 
In our perspective, these quantities can be considered as non-radial perturbations of a radially 
pulsating spacetime, which is itself described as a perturbation of a static background (see sec- 
tion l4.1l for more details). This perturbative expansion then allows us to perform a splitting of the 
perturbation in A to get 

g(l) =g (0,l) +A g(l,l) ) (55) 

where C/( 0,1 ) is a first order, gauge invariant non-radial perturbation on a static spacetime, while 
C/v 1 ' 1 ) describes the corrective effects on the non-radial oscillations brought by the radial pul- 
sations. The non-linear perturbations Q^ 1 ' 1 ^ are the quantities we have used to investigate the 
coupling between the radial and non-radial perturbations. In this chapter, we explore in more de- 
tails their structure and prove their gauge invariance with respect to a spherical static background. 
It is important to remark that the (1,1) superscript refers not only to quantities constructed from 
the g^ 1 ' 1 ) perturbations, but in general to any perturbative quantity of order Ae. This structure is 
for instance evident in equations d4.22t - d4.24t . d4.33t and in (14.6 It for some of the non-linear 
variables used in this work. Therefore, G^ ' 1 ^ and Q^ 1 ' 1 ' can be written in general as: 

g(o,i) = w (o,i) j (56) 

0(i,D = n^ + ^i^j^, (5-7) 

where the objects T^ 0,1 ) and J^' 1 ^ are linear in the (0, 1) perturbations, while ijj 1 ' ^ and H^ 1 ' 1 ) 
are respectively linear in the A and the Ae variables. 

By definition, the quantity H^ ' 1 ^ represents any of the gauge-invariant quantities given in 
the expressions d2.59t - d2.64t and d2.72b - d2.74l > for the polar perturbations, and in (I2.66t - d2.68t 
and d2.70t for the axial, which have been derived from first order non-radial perturbations of a 
static background. It is also clear that the quantities 7V 1 ' 1 ) have the same formal structure as 
the but are now constructed from second order Ae quantities. We will show later that 

these objects are not gauge-invariant at Ae order. In fact at this second perturbative order the 
gauge invariant quantities must contain some extra terms formed by the product of first order 
perturbations, as in equation d5.7t . 

Let us now consider a class of non-radial gauge transformations of first and second perturba- 
tive order where we have fixed the gauge for the linear radial perturbations. This transformation 
is given by the two expressions d5.2t and d5.4t . The linear perturbations G^ ' 1 ^ are gauge-invariant 
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by definition at order e, thus from equation (15 .2t we have that 

g(o,i) _ g(o,l) _ 7^(0,1) _ 7^(0,1) _ o . 



(5.8) 



For the non-linear perturbations we have from equation d5.7t and the fact that we have fixed 
the gauge for radial perturbations the following transformation: 



r(0,l) 



(5.9) 



This expression can be further elaborated if we note that every can be expressed 

as follows: 



(0,1) 



(5.10) 



where A and Bo- are linear operators involving differentiation with respect to the coordinates of 
M 2 and integration on S 2 . These operators act on spacetime objects and return objects with 
indices on M 2 . In equations (I5.10t and in the rest of this section we consider for simplicity only 
the metric perturbations. For the energy-momentum tensor perturbations and fluid variables the 
procedure follows along the same lines and is given in two later subsections of this chapter. Using 
the gauge transformations (I5.2t and ( 15.4b . the transformation rules for TiP-^ and ^i ' 1 ^ are given 
by the following expressions: 



£(1,1) = n^ + A 



£c s (1,0) + £ t 2 
^5(o,i)8 



(5.11) 
(5.12) 



In addition, we can notice that the term A [ £^g ] must vanish in equation (I5.11t for any vector 
field £. This result is due to the fact that the perturbative fields liP-' 1 ' and T^ 0,1 ) share the 
same functional structure, and that H^ ' 1 ^ are first order gauge-invariant perturbations. Thus, they 
satisfy equation d5.8t for any gauge transformation (I5.2t generated by a generic vector field £. As 
a result, the quantity (I5.11t becomes 

£(1,1) =H (1,D +A r 
and the gauge transformation (15 .9t reduces to the following 



£(0,1) ^ 



(1,0) 



(5.13) 



(5.14) 



The gauge transformation d5.14t is valid for Ae non-radial perturbations, when the gauge of linear 
radial perturbations is fixed. This formula will be used in the following sections in order to prove 
the gauge invariance of the C/P-' 1 ) variables. This will be achieved by showing that the right hand 



side of equation (I5.14t vanishes. 
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In order to prove the invariance of our non-linear polar and axial perturbations Q^ 1 ' 1 ' for a 
gauge transformation defined at first order by the expression fl5.2b and at second order by the 
transformation (15.4b . where the radial gauge is fixed, we proceed as follows: 



1) we expand the GSGM variables as in equation (I5.5t and determine the second or- 
der quantities C^ 1 ' 1 ). Then, we identify in QO-' 1 ) the fields TL^ l,l \ Z^ 1 ' ^ and as in 
equation (15.7b . 

2) We perform the gauge transformations of the fields TL^ 1,1 ^ and by using the rules (15.131) 
and (15.12b. 



3) Finally, we introduce the results in the gauge transformation (15.14b and prove that it is 
satisfied, namely that the right hand side vanishes. 



5.1.1 Axial perturbations 

Gauge invariant non-radial perturbations on a radially pulsating spacetime are expressed for the 
axial sector by the tensor fields (I2.66b - (l2.68b . ( 12.70b and ( 12.95b . The perturbative expansion of 
these metric and fluid variables, which has been illustrated in the previous section, leads to the 
following quantities at first order in e: 



l.(0,1) 
A 

^(0.1) 

r(0,l) 
A 

L (0,l) 



h 



(0,1) .(0,1) 



St 



(o,i) 



^0,r 



+ 2k 



(0,1) 



-(*+A) 



st m _Q /l (o.D j 



(5.15) 

(5.16) 

(5.17) 
(5.18) 



and at second order in Ae: 



£,(1,1) 

A 
*(1.1) 



h 



(1,1) 



hff+2h^v A 



r ( k 



.(1,1) 



2k, 



(i,i) 



l,t "'O.r J ' * n 
L (1,D = ( j t (l,l)_Q^(l,l)_g(l,0)^(0,l). 



-(*+A) _ (1,0)^(0,1) 



(5.19) 

(5.20) 

(5.21) 
(5.22) 



Comparing the expressions (I5.16b - (l5.18b and d5.20b - d5.22b . it is easy to identify for any second 
order perturbation the terms Tiy-^: 



1.(1.1) 
A 

$(1.1) 



_/,(!.!) +2 ^(1,D 



r ( A; 



,(l,l) 



"l ,* 



-(*+A) 



(ft 



(1,1) 



(5.23) 
(5.24) 
(5.25) 



CHAPTER 5. GAUGE INVARIANCE OF NON-LINEAR PERTURBATIONS 



67 



L (i>i) . St^-Q hP-V 
and the quantities ^^xi- 1 ' - ^J^' 1 ' 



1.(1.1) 
A 

#(1,1) 
r(l,l) 

L (1,D 



0. 



.^(1,0)^(0,1) ; 

-Q (1,0) < 1} , 
.gd,o) fe (o,i) 



(5.26) 



(5.27) 
(5.28) 
(5.29) 
(5.30) 



Now, we have to determine how the perturbative fields Tt^ 1 ' 1 ^ and Ja°^ change under a gauge 
transformation given by the expressions ( 15.131 and (15.121 . The covariant vector field £(0,1)0 that 
generates axial gauge transformations can be expanded in odd-parity tensor harmonics, 



^ 1)a = (0,0,r 2 V^S a ), 



(5.31) 



where y(°' 1 ) = V^ - 1 ) (x A ) is a scalar function on the submanifold M 2 . 

The first order metric components and M ' 1 ) change under the gauge transformation 

(15.2b as follows: 



£(0,1) 



h 



(0,1) 



+ r 2 V 



(0,1) 



A -<-' "\A 
h (0,l) + r 2y(0,l) 



(5.32) 
(5.33) 



while the master function ^o. 1 ) is gauge invariant by construction. 

The transformation (I5.13t requires the evaluation of the Lie derivatives of the radial metric 
and energy momentum tensor perturbations with respect to the generator £(o,i)a- Considering the 
expression of the radial metric in the gauge (I3.18t we have: 



p (1.0) 

X C(o,i)^«/3 



0, 



and for the energy momentum tensor 



^(1,0) 
5(o,i) AB 



St 



(1,0) 



€(0,1) ot Ab 



(1,0) 



^(0,1) "*a6 



r 2 Q (l,0) V (fi,D Sbj 

P Q m v (o,D {Sa:b + s b:a ) 



(5.34) 



(5.35) 
(5.36) 
(5.37) 



Therefore the gauge transformation ( I5.13t for the components of the metric and energy momen- 
tum tensors leads to: 



n A 



St 



2(1,1) 



A > 

6tW 



(0,1) 

A 



(5.38) 
(5.39) 
(5.40) 
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St^ = 5tf 1] + r 2Q(i,o) y (o,i) (541) 

Finally, if we introduce equations (I5.32t . (l5.33l) and (I5.38t - (l5.41t into the gauge transformation 
(I5.14t of the axial variables ( 15.191 ). (15.211) and (15.221) . we find the gauge invariance of the quantities 
L^' 1 ^ and L^ 1 ' 1 ) respectively. The gauge invariant character of the master function vj/t 1 ' 1 ) 
derives directly from its expression (I5.20t and from the fact that the radial gauge is fixed. In fact, 
the equation (I5.20t shows that it is a linear combination of the gauge invariant tensor k^' 1 ^ and 
the first order master equation ij^ ' 1 ) and the linear radial perturbation r/ 1 ' ). 

The last variable of the axial sector to study is the axial component f3 of velocity perturbation 
(I2.70t . For the first and second order gauge transformations (I5.2t and (I5.4t . the Lie derivatives of 
the background velocity and its radial perturbations vanish: 

£^a = £ ki Ju^=0, (5.42) 

where the indices (£, 1) with % = 0, 1 represent both the e linear non-radial and the non-linear Ae 
cases. Therefore, the odd-parity velocity perturbations dua' 1 ^ and 6ua^ remain unchanged and 
fiW and Z^ 1 - 1 ) are gauge invariant. 



5.1.2 Polar perturbations 

Polar gauge invariant perturbations on a radially oscillating spacetime are defined by the metric 
and fluid quantities (I2.59t - fl2.64l i. The analysis of their gauge invariance will start with the metric 
perturbations, then with the energy momentum tensor and finally with the fluid variables. The 
details of the calculations are given only for the metric variables. For the energy momentum 
tensor and the fluid variables we determine the fields TL^ 1,l \ Z^ 1 ' ^ and obtained from the 

perturbative expansion (I5.7t and the value of the Lie derivatives of the gauge transformations 
(I5.12t and (I5.13t . These quantities must be then introduced in equation ( I5.14t . where after some 
algebraic calculation the two terms of the right hand side cancel each other. 

Metric polar perturbations 

The four first order gauge invariant metric variables in the GSGM formalism are given by the 
symmetric tensor k^ B and the scalar function respectively defined in equations (I2.59E and 
(I2.60t . The perturbative expansion of these quantities provides the following expressions at first 
order in e: 

= K^»-2v A p^ 1} , (5.44) 
where the vector p^'^ i s given by the following expression: 

pT ] = - ^G\ f , (5.45) 
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and at second order in Ae: 



1.(1.1) 
K AB 

fc (l,D 



AB \Pa\B P B\A ' AB PC 

K^) -2v A p^ + 2g^ AB v B p^ 



(1,1) , (1,1) 



(1,0)C (0,1) 



,(1,0)C 



where T AB are the radial perturbations of the Christoffel symbols 



,(i,o) c 

AB 



-CD 



, (1,0) , (1,0) _ , (1,0) 
n DA\B ^ DB\A AB\D 



and the quantity p*^ is given by the following equation: 



1 



.2^(1,1) 



(1,1) (1,1) - , rv 

Pa - n A 2 \ A 



In these second order perturbations Q^ 1,l \ it is easy to identify the quantities W( 



1,1). 



(1,1) 



ft(l'l) _ fn^' 1 ' 4- n^' 1 ' 
AB ■ AB \Pa\B + Pb\A 

jfc(l.l) ; 

and the contribution brought by the first order perturbations Yln^ ^ 



(l,i) , fi,i)> 



(1,1) 
AB 



9r (l,0)C (0,1) 



fc (i,i) . ^(1,0)^^(0 



(o,i) 



(5.46) 
(5-47) 



(5.48) 



(5.49) 



(5.50) 
(5.51) 



(5.52) 
(5.53) 



Now, we consider the gauge transformations (I5.12t and (I5.13t . The generator £(o,i)a of the gauge 
transformations associated with the non-radial perturbations can be expanded in polar tensor har- 
monics, 

(5.54) 



C( ,i )a = (if 1} Y,r^Y a ). 

Under the transformation (I5.2t the components of the non-radial metric tensor change as 
follows: 



1(0,1) 
n AB 

1(0,1) 

n A 

£(0,1) 



,(0,1) 



1(0,1) 



°AB + %A\B + £ 



(0,1) 
B\A 



,(0,1) 



?(0,1) 



+ CT'+r%A 
A £(0,1) 



#(0,1) + 2C A£ 
(5(0,1) _ G (o,i) + 2^(0,1) _ 



(5.55) 
(5.56) 
(5.57) 
(5.58) 



Thus, from the definition (I5.45t and the expressions d5.55t -( l53Hl the quantity pjt , which con- 



stitutes the term JT'J ' 1 ' 1 in equations d5.52t and <I5.53I > changes as: 



(0,1) _ (0,1) p(0,l) 

Pa —Pa + ?a 



(5.59) 
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On the other hand, in the transformation ( 15.131 the Lie derivative of the linear radial metric e A 
provides the following results: 

f (1,0) _ f m,i) c. (1,0) , (i,o) p(o,i) c , (i,o) £(o,i) c\ v ( - fim 

^co,^ = t$!>FY M , (5.61) 

A ( o,i)g£' 0) = 0. (5.62) 

Therefore, the second order metric components changes with respect to the transformation ( 15.131) 
as follows: 



~h {1 ' 1] 
AB 

~h ihl) 

(5(i,D 



AB ^^^"'AB ' 



A 

G^) , 



(5.63) 
(5.64) 
(5.65) 
(5.66) 



and then the perturbation p^ 1 ^ as: 



(1,1) _ (1,1) , (1,0) ? c 

Pa —Pa + ac s ■ 



(5.67) 



Finally, we have all the elements for checking the invariance of the tensors k^'^ and fc^ 1 ' 1 ) 
under the gauge transformation (I5.14t . When we introduce in the definitions ( I5.46t and ( I5.47t the 
relative values we find: 



k 



(1,1) 
AB 



k 



(1,1) 
AB 



, £(0,1) Cl (1,0) 
^ S ABIC 



, (1,0) £(0,1) C , (1,0) £(0,1) C 
"CB SI 4 'MC S| 



"AC S| S 



B 



+ 2T 



(1,0)C£(0,1) 



AB 



c 



jfe(i,D = k dA) 



C + 2/i (1,0)AB- b |(0,1) 



(5.68) 
(5.69) 



and after some trivial calculations we get 



k A ' B ^ — k A ' B , (5.70) 
U 1 ^ = k^\ (5.71) 



thus the gauge invariance of the non-linear perturbations. 



Energy momentum tensor 

The components of the energy momentum tensor St^^ can be combined to define the seven 
quantities (I2.61t - d2.64l i. These perturbations that are gauge invariant on a spherical and time 
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dependent spacetime time can be expanded as in equation (I5.7I) . We can then write the following 
expressions: 



77(1,1) 
AB 

T (l,l) 3 
T (l,l) 2 



r + (l>l) Z _(1,1) C 

dt AB ~ t AB\CV 



T (1,1) C T (1,1) C 

tACP\ B ~ tBCP\ A 



ft(l,D 3 _ p (l,l) ( q |o + 2 g %) + + g G (l,l) > 
^' 1} -^P (1 ' 1)C " 



2 ^ U \A ' 



( y t (i,i)a_ r 2Q G (i,i) i 
and for the contribution carried by the first order perturbations Yla-^v *7tr ■ 



(5.72) 
(5.73) 

(5.74) 
(5.75) 



r(l.l) 
AS 



T (i,i) 

27(1,1) 2 



_(0,l)O # (l,0) ,(1,0) (0,1)C .(1,0) (0,1)C r a (l,0)CDjO,l) 

-tBc(g&° )CD pS ,1) +g (1 ' 0)C H°ii ) ) , (5-76) 



i(z + i; 



.Q(1,0) G (0,1) 



r 2 Q (l,0) G (0,l) 



(5.77) 

(5.78) 
(5-79) 



where Q and Qt 1 ' ) are a static scalar function and its radial perturbations, which corresponds 
for a perfect fluid to the pressure. This notation has been introduced in the GSGM formalism in 
equation ( I2.32I) . and we keep it in this section in order to distinguish the pressure from the metric 
vector pa- 

The perturbative fields H^ 1 ' 1 ^ transform accordingly to equation (I5.13t . We must then deter- 
mine the Lie derivative £^(o,i)ta^ ' where the energy momentum tensor for the radial perturba- 
tions has this block diagonal form 



One then gets 



£ 



£(0,i) l AB 



(1,0) 



£ 



5(0,1) Aa 



(1,0) 



(0,1) C ,(1,0) (1,0) £(0,1) C (1,0) £(0,1) c\ v 
l AB\C l CB $\A l AC S|b J 1 i 

(1,0)^(0,1) C j. ,.20(1,0)^(0,1)^ y 



(5.80) 

(5.81) 
(5.82) 



£ 



?(o,i) a6 



(1,0) 



Qf£' 0) l (0,1) c -l(l + l)QW>£ u > l > + 2^cQ ll ' Uj e 

+ (2r 2 Q( 1 '°)e (0 ' 1) ) Z, 



(1,0) t(0,l) c 



y 7 , 



aft 



J a6 



(5.83) 
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Therefore, the gauge transformation (I5.13t for the H^ 1 ' 1 ' terms give the following expressions: 
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(5.84) 

(5.85) 
(5.86) 
(5.87) 



The gauge invariant character of the energy momentum perturbations (I2.61t - (l2.64t at order 
Ae can be now proved by collecting the previous information. We must introduce into the gauge 



transformation (I5.14t : i) the expressions (I5.84t - (l5.87t due to the transformation of the tensor 
fields 7i^ l ' l \ and ii) the changes (I5.58t and < I5.59I > carried by the first order non-radial fields 
J"i 0,1) of equations ( B77Si -( l577g| . 



Fluid perturbations 

In the polar sector the fluid perturbations on the GSGM spacetime Mq are given by the two 



velocity perturbations (12.721) and (12.731) and energy density (I2.74t . The gauge invariant character 
of the fluid perturbation iJ^' 1 ), which is defined in equation d4.33t . will be discussed at the end 
of this section. 

The perturbative expansion ( 15.71 ) of these variables produces the following expression for the 
tensor fields H^ 1 ' 1 ': 
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(5.88) 
(5.89) 

uA 1 ) : u^-p(^ A n lA , (5.90) 
and for the part related to the nonlinear contribution of the first order perturbations 1^'°^ ji- ' 1 ^ : 
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(5.93) 



The gauge transformations d5.13t for HP-' 1 ' are now determined by the Lie derivative of the 
radial perturbations of the fluid velocity and energy density 
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(5.94) 



where Q = In p. We obtain 
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(5.95) 
(5.96) 
(5.97) 



t (l,0)|(0,l) A 

%o,d^ (1 ' 0) = i m % A . 

Therefore the gauge transformation (I5.13t for the IIS 1,1 } are 

5 (i,D = a^+u^i A -g^i A u B , (5.98) 
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(5.100) 



In addition, we can notice that in the equations d5.91t -( l5T9"3l the first order non-radial terms 
J^ ' 1 ^ are given by the metric quantities h^B jPa an d tne fluid velocity 5u^'^ ■ The behaviour 



of h AB ' and p A ' ' under the gauge transformation d5.2t is expressed by the equations (15.551 
and d5.59t . while the velocity perturbation changes as follows: 



Su 
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(5.101) 



We can now bring the transformations d5.98b - (l5. lQOt . d5.55b . d5.59t and d5.101t into the gauge 
transformation (15.14l >. then find the gauge invariance of these non-linear quantities. 

Finally, we discuss the fluid perturbation (14.321 . Its perturbative expansion provides the 
following expressions at first and second order: 



#(0,1) = n (o,i) = yP_ jo,i) 

p + p 



(5.102) 
(5.103) 



where we identify 
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= J ' 1 ). (5.106) 

Hence, the gauge-invariant character of H^ ' 1 ^ and H^ 1 ' 1 ' for a gauge transformation ( I5.14t with 
a fixed radial gauge follows from the gauge invariance of a/ ' 1 ) and a/ 1 ' 1 ', previously proved (see 
equation (I5.93t ). 



Chapter 6 

Numerical Simulations 



This chapter is dedicated to the description of the numerical code that simulates the non-linear 
dynamics arising from the coupling between radial and non-radial perturbations. In particular, 
this numerical analysis is focused on the axial non-radial oscillations, as for the polar non-radial 
perturbations the implementation of the code is still under way. The linear and non-linear per- 
turbative equations form a hierarchical boundary initial value problem, where initial values for 
the linear radial and non-radial perturbations can be independently set up. The two independent 
initial configurations that we can investigate are given by: i) a radially pulsating and differen- 
tially rotating star, and ii) the scattering of an axial gravitational wave on a radially pulsating star. 
The former configuration is the more interesting. In fact, at first perturbative order the axial non- 
radial perturbations of the fluid quantities do not have any dynamical properties. The only matter 
perturbation is given by the axial velocity that describes a stationary differential rotation. As a 
result, the star is not a source of gravitational radiation. This aspect changes radically at second 
order, when the radial oscillations couple with the differential rotation. Now, the presence of this 
stationary axial velocity and the related frame dragging allows the radial perturbations to exhibit 
their pulsating character and drive the oscillations in the source terms of the second order axial 
master equation. This non-linear oscillating dynamics then produces an axial gravitational signal 
that as we will see in section 16. 5 . 3 1 precisely mirrors at coupling order the spectral properties of 
the radial perturbations. In addition, it is worthwhile to remark that this is a first order effect with 
respect to the axial velocity perturbation and that it is strictly related to its differential character. 
In fact, it is well known that for a uniform rotating star Q = const the quasi-radial modes appear 
at second perturbative order in 17, when the deformation of the star is taken into account by the 
perturbative analysis. In section 16.11 we introduce the structure of the numerical code, focusing 
on its hierarchy, implementation of the numerical grids and analysis of characteristic curves. In 
addition, we describe also the introduction of the "tortoise fluid coordinate", which provides a 
more accurate description of the radial pulsations near the stellar surface. The background solu- 
tion is illustrated in section 16.21 while the radial and axial non-radial perturbations are described 
respectively in sections 16.31 and section 16.41 In these sections, we write the numerical methods 
for simulating their evolution and discuss in detail the setting up of the initial configurations. Fi- 
nally, section 16.51 is dedicated to the non-linear simulations for the description of the coupling 
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between the radial and axial non-radial perturbations, where we provide the numerical methods 
and discuss the results. 

6.1 Numerical framework 

The numerical code for the evolution of non-linear oscillations in the time domain reflects the 
hierarchy structure of the non-linear perturbative theory. Starting from a solution of the TOV 
equations that describes the equilibrium configuration of a nonrotating spherical star, we must 
solve independently the two classes of first order perturbations, i.e. the radial and the axial non- 
radial for an arbitrary harmonic index I. These first order values are necessary at any time step for 
updating the sources that drive the non-linear oscillations. Eventually, the perturbative equations 
for the second order perturbations can be integrated (see figure l6~Tt . 

The separation in the perturbative fields of the angular dependence from the other two coor- 
dinate, which derives from the 2+2 splitting of the four dimensional spherical spacetime and the 
existence of a tensor harmonics basis for the 2-sphere, leads to a 1+1 dimensional problem. One 
dimension is given by the spatial coordinate r and the other one by the coordinate time t. 

The numerical simulations in this code are based on the "finite differencing method", where 
derivatives are replaced with finite difference approximations (see Appendix [FJ. As a result, the 
differential equations become a set of algebraic equations that can be solved with standard numer- 
ical methods for partial and ordinary differential equations [91 ]. In particular, we adopt explicit 
numerical schemes with second order approximations for the spatial derivatives and first or sec- 
ond order approximations for the time derivatives. The use of a first or second order numerical 
algorithm in time will depend on the particular properties of the perturbative equations. 

The numerical stability and dissipation of the simulations will be monitored with the L2- 
norms of the numerical solutions and its accuracy with a convergence test, see Appendices (IG.31 
and (IG.21 respectively. 

Fluid tortoise coordinate transformation 

The radial perturbations can manifest a reduction in accuracy when the speed of sound tends to 
vanish, near the stellar surface, where the characteristic curves of the sound wave equations (13.291 
and ( 13.441 and of the radial perturbative system of equations ( 13.241 and (13.251 . which are given 
by the following expression: 

£i = r ± v s t where v s = c s e* _A , (6.1) 

lose their propagation character (see figure lo3l ). In practice, the numerical simulations carried 
out with a second order scheme show that the convergence rate of the radial perturbative solutions 
falls to one immediately after the wave packet touches the surface. This behaviour has been well 
described by Sperhake in his PhD thesis [ 101 ], here we report only the important aspects neces- 
sary for our work. A method for solving this problem is given by a refinement of the numerical 
grid. However, since the accuracy issues arises near the surface it is not necessary to increase the 



CHAPTER 6. NUMERICAL SIMULATIONS 



77 



Background 
TOV 



Initial Values 
Radial Pulsations 



Radial Pulsations 
Evolution 



Initial Values 
Non-Radial Oscillations 



Non-Radial Oscillations 
Evolution 



Sources 



Non-linear Coupling 



Figure 6.1: Code hierarchy for the time evolution of the non-linear axial perturbations arising from the 
coupling between the radial and non-radial oscillations. The initial configurations of the first order radial 
and non-radial perturbations are independent. 



resolution homogeneously on the whole grid. It would be more appropriate to perform a coordi- 
nate transformation that, close to the surface, simulates a high refined grid for the r coordinate. 
This characteristics are satisfied by the "tortoise fluid coordinate", which have been introduced by 
Ruoff for the analysis of stellar non-radial perturbations where the matter is described by realistic 
equation of state l95ll . 

The new spatial variable x in the "tortoise fluid coordinate" is defined as follows: 

dr = c s dx (6.2) 

where c s is the speed of sound. The decreasing character of this velocity (figure I6.4t and the 
definition (I6.2t imply that an evenly spaced grid with respect to the new coordinate x is able to 
simulate a grid for r which will be more and more refined toward the surface. 

This aspect is particularly important at second perturbative order, where the perturbative equa- 
tions contains source terms in the interior spacetime. A coarse resolution in the low density re- 
gions near the surface could produce some spurious oscillations in the radial pulsations, which 
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could be propagated into the exterior spacetime through the junction conditions and then affect 
the gravitational signal. 

Therefore, we introduce the coordinate transformation d6.2t in the part of the code dedicated 
to the radial pulsations and consequently to the TOV equations that describe the hydrostatic equi- 
librium. In practice, we must replace the spatial derivatives with respect to the coordinate r with 
the following expression: 

d r = —d x . (6.3) 

c s 

The new equations for the radial perturbations and for the background are given in sections 16.21 
and [O] respectively. From the radial perturbative equations ( 16.301) and ( I6.31t we can notice that 
with the tortoise coordinate the velocity of the characteristic curves 

£i = x ± v s t where v s = (6.4) 

does not vanish on the surface. 

6.1.1 Numerical grids 

The fluid tortoise coordinate transformation concerns the radial perturbative fields which are 
present in the interior spacetime only. On the other hand, the linear and non-linear axial per- 
turbations are well described by the Schwarzshild coordinate r on the entire spacetime. As a 
result these two perturbative families are defined on two distinct integration domains that must 
be maintained separate in the construction of the numerical code. These domains are both 1+1 
dimensional, where one dimension is given by the time coordinate t and the other by the spatial 
coordinate x for the radial perturbations and r for the non-radial perturbations. 

The 2-dimensional and continuous evolution domain V r C M 2 for the non-radial perturba- 
tions is discretized along the two dimensions (i, r) with an evenly spaced mesh that we call the 
r-grid: 

rj = n+jAr with j = 0, l,..,J r , (6.5) 
t n = ti + nAt with n = 0,l,..,N r , (6.6) 

where the quantities Ar and At are the constant spatial and time increments, and J r and N r 
denote the number of grid points. The background and perturbative fields of the linear T^ ' 1 ) and 
non-linear Ti 1 ' 1 ) axial oscillations are then approximated by a set of discrete quantities evaluated 
on the points of the numerical grid 

(T^y. ee T(W {tnirj) t ( rW )); ^ T W) (t n , rj ) , (6.7) 

where the upper index n denotes the time level and the lower index j the spatial mesh point. We 
have shown in equation (I6.7t only the perturbative fields, for background quantities the discrete 
approximation is obviously similar. 
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The 2-dimensional and evolution domain V x C M 2 for the radial perturbations is instead dis- 
cretized along the two dimensions (t, x) with an evenly spaced mesh, the x-grid: 

Xj = xi + jAx with j = 0, 1, .., J x , (6.8) 
t n = nAt with n = 0, 1, ..,N X , (6.9) 

where Ax and At are now the constant increments for the tortoise fluid coordinate and the time, 
which is in general different from the representation given in the r-grid. As before, the integers 
J x and N x denote the grid dimensions. The radial perturbative fields 7~( 1,0 ) are then discretized 
on the x-grid as follows: 

(T^) n ^T^(t n ,x,). (6.10) 

From the definition ( I6.2t of the tortoise coordinate, we can see that the radial simulations carried 
out on the x-grid, can also be mapped on this new representation of the coordinate r: 

fj = fi+jAr where Af = c s Ax and j€N, (6.11) 

which for a polytropic equation of state has an increasing resolution forwards the surface. 

Interpolation 

The implementation of the two grids in the numerical code is shown in figure 16.21 The TOV 
solutions for the equilibrium configuration are discretized on both the x- and r-grids. The simu- 
lations of linear radial perturbations are carried out on the x-grid while the linear non-radial on 
the r-grid. Therefore, in order to provide the source terms quantities evaluated at the same spatial 
mesh points, the radial perturbations are interpolated on the evenly spaced r-grid. In particular, 
this procedure is constructed between the radial quantities determined in the non-homegeous rep- 
resentation f d6.11t of the spatial coordinate r and the r-grid. With the updated values of the 
source terms the non-linear simulations can be then carried out. More precisely, due to the im- 
plementation of explicit numerical schemes the evaluation of the non-linear perturbations at the 
n + 1 time slice will rely on the source determined at the previous time slice n. This property 
can be seen directly from the discretization schemes given in the following sections. In addition, 
the interpolation displays more accurate results when the x-grid dimension is twice the r-grid, 
namely J x = 2J r . 

In the next section, we discuss the Courant-Friedrichs-Levy condition and show how to set up 
numerical simulations on the two grids adopting the same Eulerian time. 

6.1.2 Characteristic curves and Courant-Friedrichs-Levy condition 

A discrete representation of the spacetime must satisfy the Courant-Friedrichs-Levy (CFL) con- 
dition which is a necessary but not sufficient condition for the stability of the code. In order to 
well describe the evolution of a physical system and its causal structure, the numerical domain 
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Figure 6.2: Integration scheme of the numerical code which simulates the coupling between the 
radial and non-radial perturbations by using two different one-dimensional grids and an interpo- 
lation procedure. 



of dependence must include the physical domain of dependence. This implies that the physical 
velocity v of the system has to be lower than the numerical velocity, i.e. 

Ar 

v<—. (6.12) 

This equation can be also used to determine the maximum time increment allowed for a given 
velocity and a spatial resolution, 

Ar 

At max = . (6.13) 

V 

In this code, we have set up the time and space grid increments in order to: i) satisfy the CFL 



condition (I6.12t and ii) have the same discrete representation of the Eulerian time coordinate, 
namely At = At. This latter requirement allows us to couple, in the source terms, the first order 
perturbations which are evaluated at the same instants of time. 

The physical velocity v gw of the axial gravitational waves can be determined by the charac- 
teristic curves of the master equations ( 13.751 and d4.62t . For the radial perturbations, which are 
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Figure 6.3: Characteristic curve velocity profiles v gw for the master axial equations J3.75i and ( I4.62t . 
and v s for the radial hyperbolic system of equations (l6.30t - d6.3U . The radial propagation velocity is plotted 
in x-coordinate. The stellar radius is at R s = 8.862 km in the r-coordinate and is at R x s = 25.80 km in 
the ^-coordinate. 



evolved on the x-grid, the sound velocity v s is given by the system of equations d6.30t - (l6.31t . 
They read 

v s = e*" A for ie[0,i?], (6.14) 
v gw = e*" A for re[0,oo), (6.15) 

where the radius of the star with respect to the x coordinate has been defined by the quan- 
tity R x s . The two numerical velocities are formally the same, but v s is defined on the x-grid and 
is present only in the interior spacetime, while v gw is the velocity of the gravitational wave in the 
Schwarzshild coordinate and is defined in the whole spacetime. The profile of these velocities for 
our equilibrium stellar model is plotted in figure l631 where the stellar radius is at R s = 8.862 km 
in the r-grid and is mapped by the tortoise coordinate transformation at R x s = 25.80 km in the 
x-grid, (see section l6~2l . If we define the numerical velocities of the r- and x-grids respectively 
as follows: 

Vr = -IT- : V T = , (6.16) 

At ' At 

the CFL condition (I6.12t implies that 

v s < v x v gw < v r . (6.17) 
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Therefore, in order to have the same Eulerian description of the time coordinate in the two grids 
the following relation between the numerical velocities must be imposed: 



We set up in this code the two numerical r- and x-grids for the interior spacetime, where 
the x-grid dimension is twice the r-grid, J x = 2J r . The coarse dimension for the r-grid is 
J r = 200, which for the stellar model adopted in this thesis leads to the spatial increment 
Ar = 0.04431 km. This choice leads to an x-grid with a dimension J x = 400 and spatial 
increment Ax = 0.0646 km. The CFL condition will certainly be satisfied if we fix the numeri- 
cal velocity of the r-grid as v r = 0.99. In fact, the corresponding value for the x-grid velocity 
is v x = 0.67905, which has been determined from equation (I6.18t . This value is higher than the 
physical velocity v s shown in figure 16.31 for the stellar model adopted in this work. The same 
properties are valid also for higher resolution for the x-grid as long as the relation J x = 2 J r is 
maintained between the two meshes. 

In the exterior spacetime only the r-grid is present which conserves the spatial and time steps 
of the internal mesh. 

6.2 Background 

The background spacetime is represented by a perfect fluid, spherically symmetric relativistic 
star in hydrostatic equilibrium. The TOV equations d3.5l >-( l3~Hl and an equation of state, which 
describes the properties of the stellar matter, form a closed system of equations, which can be 
integrated by specifying the central density of the star. As explained in section |3~T1 in this work 
we will consider a poly tropic equation of state ( 13.9ft . 



During the numerical integration of the TOV equations, we have found a slight improvement in 
the rate of convergence by adopting the metric function A as a new independent variable instead 
of the mass function M. Furthermore, the presence of two numerical meshes in the code requires 
an integration of the TOV equations in both the r- and x-grids. 

The equilibrium configuration is then the solution of the following system of equations: 



p = k p T . 



(6.19) 




(6.20) 



(6.21) 



D<£> 



Tip 



(6.22) 



P + P ' 



Dr 



1 



(6.23) 
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p = Kp r , (6.24) 

where we have introduce the differential operator D, which acts on a generic scalar function / as 
follows: 

~ \f r for the r-grid , 

D/ = < (6.25) 
I c s f,x for the x-grid . 



This definition allows us to write the same expressions for the TOV equations ( I6.20I) - (I6.23I) for 
the "tortoise fluid coordinate" (I6.2t as well as for the (area) coordinate r. By adopting this new 
system of equations the mass function M can be determined in terms of the metric function A 
with its definition (I3.4I) . In case of an integration carried out in the "fluid tortoise coordinate" 
frame, the equation ( 16.231 ) shows that the radial coordinate r is a scalar field which is an unknown 
of the system as well as the other metric and fluid variables. 

For specific values of the polytropic parameters K and V in the EOS J6.24I) . the numerical 
integration of equations ( I6.20I )-( |6~2"4T ) can proceed from the stellar origin r = outward as fol- 
lows 1 77]: 

1) Specification of the central mass-energy density p c and consequently the corresponding 



value of pressure p c determined trough the EOS (I6.24I) . 

2) Imposition of the boundary conditions at the origin of the coordinates r = for the metric 
variables <E> and A: 

$(0) = S C , A(0) = 0. (6.26) 

The second condition in equation ( 16.261 ) can be derived from the definition ( 13.4b . and from 
the behaviour M ~ 0(r 3 ) given by equation d3,8t . On the other hand, the choice of $ c is 
completely arbitrary. Due to the linearity of the ODE ( I6.22I ). its value can be later rescaled 
to the correct one in order to satisfy the matching condition on the stellar surface. 

4) Integration of the equations ( I6.20l i-( l6~2"4l with a standard numerical method for ODEs [91 ]. 
The integration starts from the origin onward and ends when the pressure vanishes p(r) = 
0. This position of null pressure identifies the stellar surface and then the stellar radius R s . 
Thus the value of the mass function in this point M = M(R S ) is the total gravitational 
mass of the star. 

5) Eventually, the function $ can be rescaled with a constant value in order to match the 
Schwarzschild solution on the surface 

$(Rs) = -A(Rs) = \ In (l - \ . (6.27) 



We consider a stellar model which is described by the polytropic EOS (16.241) with the following 
parameters 

K = 100 km 2 T = 2 . (6.28) 

This choice allows us to determine a star with a mass and radius similar to those obtained with 
realistic equation of state of average stiffness. For a supranuclear central density given by p c = 
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Figure 6.4: For a poly tropic non-rotating star with indices r = 2, k = 100 km 2 and with a central 
density p c = 3 x 10 15 gcmT 3 ', we plot the spatial profile of {clockwise from top left): the mass energy 
density p, mass function M, speed of sound c s and the pressure p respectively. The plots are against the 
r-coordinate (solid line) and x-coordinate (dashed line). 

3 x 10 15 gcm~ 3 , the integration procedure carried out with a fourth order Runge-Kutta (RK4) 
method [91] provides a star with a radius of R s = 8.862 km and a mass of M = 1.869 km = 
1.266 Mq. When we use the x-fluid tortoise coordinate, the stellar radius is mapped at R x = 
25.840 km. In figure 16.51 one can notice the relation between the f and x representation of 
the spatial coordinate r and the higher point density near the stellar surface. The equilibrium 
configuration for all the fluid and metric quantities are plotted in figures (I6.4t and <6.5b . In the 
exterior the Schwarzschild solution is described only by the metric functions $ and A that are 
written in terms of the total gravitational mass of the star as follows: 

$(r) = -A(r) = iln(\-^ , (6.29) 

where their values on the stellar surface are given by the junction condition (I6.27t . 

Finally, we have determined the rate of convergence of the TOV solutions with the method 
described in appendix IG.2l bv using three grids of 200,400, and 800 points. The results are shown 
in table 16. II where a r and a x are the convergence factor for solutions determined in the r- and x- 
coordinates respectively. The fourth order convergence is found in accordance with the numerical 
scheme used (RK4). 
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Figure 6.5: For the same stellar model adopted in figure loT4l the metric quantities <£> and A are shown 
in the two upper plots, where the solutions determined in the r-grid are denoted with the solid lines 
while the dashed lines are relative to the x-grid solutions. In the lower figures, the relation between 
the Schwarzschild /--coordinate and the fluid tortoise x-coordinate is shown on the left, while the sound 
velocity v s is on the right. 



Convergence 
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A 




4.017 


4.007 


4.007 


4.007 


4.046 


4.031 


0~ r 


4.055 


4.000 


3.997 


3.995 


3.987 


4.055 



Table 6. 1 : Convergence rate of the TOV solutions for the stellar model considered in this thesis, where a x 
and ay are calculated for the x- and r-grids respectively. The convergence factor has been calculated for a 
set of three grids of dimension (200,400,800). The results confirm the fourth order convergence expected 
by the accuracy of the RK4 method. 

6.3 Linear radial pulsations 

The time-evolution of the linear radial pulsations of a static star can be studied with the system 
of equations (13.24t - (13.28l > for the four variables S^- ,0 \ rj^ 1 ' ' , j^' \ H^- ,0 \ which has been pre- 
sented in section (I3.3.1t . The presence of hyperbolic and elliptic partial differential equations in 
this system allows us to choose two integration approaches: 

1) a purely hyperbolic formulation (PHF), where we solve the three evolution equations (I3.24t . 
(I3.25t and (13.261 1 and we monitor numerical errors in time by looking at the Hamiltonian 
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2) A hyperbolic-elliptic system of equations (HEF), where we integrate Eqs. (13.241) and (I3.25t 
and then solve the Hamiltonian constraint (I3.28t for the metric variables S^ 1,0 ). In this case 
the constraint will be satisfied by construction. 

As shown in figure l6~2l the radial perturbations are first integrated on the x-grid. These solutions 
are then interpolated on the r-grid at any time slice. In terms of the fluid tortoise coordinate x de- 
fined in equation d6.2t the system of perturbative equations (I3.24I) - (I3.27I) becomes the following: 
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Xt 
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-c s e ^ c s 
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4irp r + + - e 
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(6.30) 



H(i,o) _ 4vr { p + p) r e *+A H (i,o) _ ( 47 rpr 2 + ± ) e *+ A s(*.°) , (6.31) 
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-Mp + p) e * +A 7 (1,0) 
c s Att(p + p)r 



„2A 



(6.32) 
(6.33) 



and the Hamiltonian constraint (I3.28t is given by: 

<?(i,o) = - c ( 87T p r -l + 2 ™\ e 2A 5 (i,o) + 8n P+l e 2A #(1,0) (6 34) 



6.3.1 Numerical algorithm 

The two evolution equations (I6.30t and (I6.31t for the enthalpy flA 1 ' ) and the radial velocity 
^(1,0) perturbations can be solved with various PDE numerical algorithms [91 ]. We have used a 
two-step McCormack algorithm (see appendix [GJ, which implements a predictor and corrector 
step and provides second order approximation in space and time. For the predictor step the finite 
difference approximation of equations (I3.24t and (I3.25t reads: 



H 



n+l 



At 



(6.35) 



#n 



+At (6 3 ). 



era 1 era 



(6.36) 



The values found, which have been denoted with a tilde, are now used in the corrector step as 
follows: 
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3,3 Ax 2 



-At (b 



(6.37) 
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(6.38) 
(6.39) 



The coefficients (b\)j, 
ground quantities: 

h = - 

b 2 = 
b 3 = 



and (63 )j are the discrete approximations in rj of the following back- 
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-4tt (p + p) re 
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(6.40) 
(6.41) 
(6.42) 



The metric variable S^ 1 ' ) can be updated at every time-step with one of the two equa- 
tions (I6.32t and d6.34t . In a purely hyperbolic formulation (PHF) the variable is evolved by 
the equation (16.32b . which can be solved with an up-wind method [91 1 where we have introduced 
a second order approximation in space. The numerical algorithm is then given by 



S? + At (64) 



n _i_ n 

Ij +7j+i 



(6.43) 



where the coefficient 64 is 



b A = -8vr (p + p) e* +A . (6.44) 
Alternatively, in the hyperbolic-elliptic formulation (HEF) discussed above we can integrate the 



Hamiltonian constraint ( I6.34t as an ODE at any time step. The equation ( I6.34t can be discretized 
by using a second order finite approximations in space and then written as a tridiagonal linear 
system. We can then use a standard LU decomposition and a "tridiagonal subroutine" |91] for 
getting the value of S^'°\ The components of the LU decomposition of equation (16.341) are given 
by the following expressions: 



a 3^~l + fyS™ + c j&j+l — fj 

where the coefficients a,b,c, and / are 
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(6.45) 



2 2m 



-c s 8npr h 



„2A 



-2A j 



2_ A_X ' 



(6.46) 
(6.47) 

(6.48) 
(6.49) 



CHAPTER 6. NUMERICAL SIMULATIONS 



88 



The other metric perturbation 7/ 1 ' ) can be solved by integrating equation ( I6.33t . At any time 
slice, this equation is a two point boundary value problem that must satisfy the condition (I3.33t at 
the origin and vanishes on the surface. We use a shooting method, where the integration is carried 
out from the origin to the surface by specifying an initial guess for r/ 1,0 ) at origin. The algorithm 
will shoot the solution up to the surface and control whether the surface condition is satisfied. If 
not, the routine will correct the value of the variable at the origin and repeats the operation until 
the surface condition is fulfilled. 



6.3.2 Boundary conditions 



The condition at the origin for the radial perturbations are given by equations d3.32l i- J3.35l l. We 
implement a grid where the first point of the spatial coordinate r\ is not at the origin tq = 
but at r\ = ro + Ar, where Ar is the spatial grid step. The boundary condition at the origin of 



the perturbative variables is then implemented at n by using the behaviours (I3.32t - (l3.35t . which 
have been obtained with a Taylor expansion. This procedure is particularly useful in the axial 
sector ll95l for avoiding numerical instability which can be generated by the presence of the term 
1(1 + l)r~ 2 in the potential of the axial master equation [95 ]. This method for the radial quantities 
7 ( 1 '°) and gives 

7? = -72, Si n = -S 2 n , (6.50) 

T2 r 2 



where we have assumed that the behaviours (I3.32t -( l3.35t are valid for both the first and sec- 
ond grid points r\ and r%. For the enthalpy i2V>°) and the other metric perturbation 7/ 1 ' ) we 
implemented the following conditions: 

at-***:**, tf-^Sf*, (6.M) 



which have been determined by spatial differentiation of equations (13.341 and (I3.33I) . and using 
the second order finite one-sided approximation ( IF.7I ) for the first order derivative. This condition 
for r/ 1 ' ) can be also derived by the constraint (I6.33t . 

At the surface the vanishing of the Lagrangian perturbation of the pressure leads to equa- 
tion (I3.42t for the variable The stellar matter in our model is described by a polytropic 
EOS, where pressure, mass energy density and speed of sound vanish on the surface. Therefore, 
the condition (I3.42t is certainly satisfied if the velocity 7( 1,0 ) and its spatial derivative 7,^'°^ are 
finite. According to the fluid tortoise transformation ( 16.31 ) the spatial derivative in the tortoise 
coordinate is related to the derivative in r by the following expression: 

7&°> = csl^ • (6.52) 

Hence, on the stellar surface we can impose the vanishing of 7, x which also ensures the finite- 
ness of 7C 1 ' ) and 7,r 1951 . The second order finite one-sided approximation (IF. 81) of the first 
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Figure 6.6: Eigenfunctions of the radial perturbation ^ The upper panel displays the eigenfunctions 
of the fundamental mode and the first three overtones, while the lower panel from the fourth to the seventh 
overtones. 



order spatial derivative leads to the following expression: 

475.-1 -7?,- 2 

' Jx O 



(6.53) 



The other values for the enthalpy and the metric variable S^ 1 ' ) can be directly determined by 
the perturbative equations d6.30t and (16.34b . while for r/ 1 ' ) we can easily implement its trivial 
condition, i.e. rj^ 1 ' ' = 0. The resulting expressions are the following: 



<-> 7 

" X 
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Jx 



on 



0. 



(6.54) 
(6.55) 
(6.56) 



where we have used the finite difference approximation (IF. 81) also for the condition (I6.54I) . The 
coefficient 65 is so defined: 



4npr + 



771 



(6.57) 
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6.3.3 Initial configuration for radial pulsations 

The time-domain integration of the radial perturbative equations d6.301 -( l633l requires the choice 
of initial values for the variables 7(°> 1 ), iJv ' 1 ), S^ - 1 ) and r/ ' 1 ) on aCauchy surface. However, the 
presence of two constraint equations implies that we can independently specify two perturbations. 
In addition, the radial perturbations must satisfy the boundary conditions on the initial slice as well 
as along the time evolution. These conditions are certainly satisfied if we specify a profile only 
for the radial velocity perturbation 7( 1,0 ) and set to zero the enthalpy perturbation H^ ' 1 ^, as the 
constraint equations (16.331 and (16.341 and boundary conditions imply the vanishing also of the two 
metric perturbations S^ ' 1 ) and r/ ' 1 ). In the case of a single eigenfunction, this corresponds to a 
choice of the origin of time. Indeed, as we can see from the radial perturbative equations (16.301 - 



(I6.33l l. we can consistently choose a normal mode oscillation with eigenfunction u) n to have the 
form: 

7 d,o) = 7 (i,o) (r)cosWnt) #(i,o) = H W\r)sinu n t, 

(6.58) 

5(1,0) = s£>°Xr)wi.w n t, v (m =??i 1,0) (r)smw n i, 

thus at t = only 7c 1 ' ) is nonvanishing. With this choice the Hamiltonian constraint (I6.34t is 
initially satisfied by construction. This choice is also consistent in the case of an initial Gaussian 
pulse, which can be considered as a particular linear conbination of the eigenmodes (16.581 . 

The time evolution of radial pulsations of a spherically symmetric perfect fluid star is com- 
pletely determined by a complete set of radial eigenfunctions. In this work, we consider two 
different configurations for the dynamical evolution, respectively described by: i) the presence 
of a broad range of normal mode frequencies, ii) the excitation of a few selected radial modes. 
The first configuration can be accomplished by imposing for 7c 1 ' ) an initial Gaussian pulse so 
that many radial modes are excited at the same time. This kind of initial data is mainly used to 
test and calibrate the reliability of the code. On the other hand, the second configuration can be 
determined by imposing as initial data the eigenfunctions of yQ>°) in order to excite the desired 
associated eigenfrequencies. 

The eigenvalue problem for the radial perturbation 'yO-fii has been described in section rj.3.31 
However, since the simulations for the radial perturbations are carried out on the x-grid we have 
to introduce the tortoise fluid coordinate also in the Sturm-Liouville problem (13.511 . The system 



of two ordinary differential equations d3.521 - d3.531 transforms as follows: 

yg' 0) = csj'°\ (6-59) 

z%® = -c s (u 2 W + Q)y^, (6.60) 

and the functions W, P, Q are given by: 



r 2 W = (p + p) e 3A+<I \ 



(6.61) 



CHAPTER 6. NUMERICAL SIMULATIONS 



91 



Normal mode 


Frequecy domain 
[kHz] 


Time domain 
[kHz] 


Kokkotas and Ruoff 
[kHz] 


Relative error 


F 


2.138 


2.145 


2.141 


0.10 % 


HI 


6.862 


6.867 


6.871 


0.13 % 


H2 


10.302 


10.299 


10.319 


0.16 % 


H3 


13.545 


13.590 




0.33 % 


H4 


16.706 


16.737 




0.18 % 


H5 


19.823 


19.813 




0.05 % 


H6 


22.914 


22.889 




0.11 % 


H7 


25.986 


25.964 




0.08 % 



Table 6.2: The table shows the eigenfrequencies of the first eight normal modes of the radial perturbations, 
which have been determined with the Sturm-Liouville problem (second column) and with an FFT of the 
time evolution (third column). The first three normal modes obtained in frequency domain have been 
compared with those published by Kokkotas and Ruoff [ 59 1 (fourth column), and the relative error is given 
in the first three rows of the last column. The remaining rows display the relative errors between the 
frequency determined in the time and frequency domains. 
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(6.62) 
(6.63) 



Now, the boundary conditions at the origin and surface are: 



- Z Q 

yo = c s — 



(p + p)c s e-*yW 



r=R x 



0. 



(6.64) 



The method used to integrate the radial eigenvalue system of equations d6.59l >- <l6.60l > is the "re- 
laxation method" |9*T1 . The first eight eigenfunctions for the radial velocity ^t 1 - ) are plotted 
in figure 16.61 with respect to the r coordinate and after a normalization with the norm of their 
maximum value. These profiles have the characteristic node numbers expected by the theory, 
i.e. absence of nodes for the F-mode, and a node number equal to the order of the overtone. The 
associated eigenfrequencies are written in table !6.2l The first three have been compared with pub- 
lished values obtained by Kokkotas and Ruoff [ 59 1 for our stellar model. The results are accurate 
to better than 0.2 percent (see table I6.2t . The rate of convergence for the eigenfrequencies and 
eigenfunctions is of second order, as expected by the accuracy of the numerical method. 



6.3.4 Simulations of radial perturbations 

The analysis of the radial oscillation part of the code starts with a pulsating configuration de- 
scribed by selected radial modes. The desired oscillation frequency is excited by introducing the 
initial condition (16.58b . where at t = the eigenfunction of the velocity perturbation -y^ 1 ' ) is 
given by 

7£' 0) = A(^ 0) (6-65) 

where we have introduced a constant factor A^ 1 ' ^ to control the amplitude of the oscillations 
and jn'°^ is one of the normalized eigenfunctions plotted in figure 16.61 The first simulation is 
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Figure 6.7: Time evolution of the four radial perturbations, where the oscillations have been excited with 
the F-mode eigenfunction of the variable -y^ 1,0 ) . The quantities have been averaged in the interior spacetime 
by using the definition (I6.66> . 



carried out for the fundamental mode with an initial amplitude of A^ 1 ' ' = 0.01 and for an evolu- 
tion time of 4 ms. The x-grid has dimension J x = 400, while the time step is chosen as explained 
in section 16.1.21 in order to satisfy the CFL condition and have the same Eulerian time discrete 
representation of the r-grid. The evolution of the four radial perturbations j0->°\ H^'°\ S^ 1 ' ) 
and j/ 1 ' ) is shown in figure l6~71 where we have plotted the following average values determined 
at any time step: 

1 

<f>= ir fdr, (6.66) 

Ms Jo 

in order to have global information about the oscillation properties of the variables studied. The 
function / in equation (I6.66t obviously represents one of the four radial perturbations cited above. 
The results show the typical periodic character of the adiabatic radial pulsations. This monochro- 
matic character is also confirmed by their spectra (figure l6~T2l which have been determined with 
a Fast Fourier Transformation (FFT) of the simulations. 

In order to test the stability and possible numerical dissipative effects, we perform a longer 
simulation of 30 ms and we monitor the Hamiltonian constraint and determine the L2 norms of 
the variables under consideration. In figure l6~8l we show the oscillation amplitude on a logarith- 
mic scale of the four radial perturbations, i.e jt 1 ' ) , ij( 1 >°) , gt 1 ' ) ( r}^'°\ and the numerical errors 
due to the violation of the Hamiltonian constraint for a hyperbolic-elliptic formulation (HEF) 
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Figure 6.8: The fundamental mode oscillations of the four radial perturbations are compared with the 
numerical errors due to the violation of the Hamiltonian constraint. The quantities are plotted in semilog- 
arithmic scale after having performed the spatial average as defined in equation (16.661 . In the upper and 
lower panels are shown the results for a HEF and PHF formulation respectively. In both figures the curves 
relative to the Hamiltonian constraints are represented in {solid lines) while the radial perturbations in 
dashed and point-dashed line. For the details of the results, see the discussion in section lo". 3. 41 

{upper panel) and a purely hyperbolic formulation (PHF) {lower panel). In the HEF, the Hamil- 
tonian constraint {solid lines) remained bounded and is three orders of magnitude lower than the 
amplitude of the radial perturbations j^ 1 ' ^ and H^'°\ which are denoted with a dashed and 
point-dashed line respectively. In addition, when we perform the average (I6.66t of the Hamilto- 
nian constraint by neglecting the first two grid points 7*1 and r-i, we find an appreciable reduction 
of the numerical errors. As shown by the lowest solid line the numerical errors are less than 10~ 18 . 
This accuracy is expected as the Hamiltonian constraint in the HEF is solved by updating one of 
the radial variables. However, the results show that the implementation of boundary conditions 
at the centre of the star introduce some numerical errors which is anyway three orders of magni- 
tude lower than the radial physical oscillations and is not propagated along the star. We perform 
a similar analysis for the PHF (lower panel of figure I6.8t by comparing the two metric radial 
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Figure 6.9: Norms of the radial perturbations ^ 1 '° s > and H^ 1 ' ^ for a 30 ms simulation where only the 
F-mode has been excited. 



perturbations S^ 1 ' ' and rf 1, °' {dashed and point-dashed line respectively) with the numerical os- 
cillations due to the violations of the Hamiltonian constraints {solid line). We can conclude that 
even for this case the constraint is well satisfied. 

The stability of the radial simulations is confirmed also by the analysis of the L2 norms. For a 
30 ms evolution excited by the eigenfunction of 7( 1 >°) associated with the F-mode, the L2-norms 
for i( 1,0 ) and H^ 1 ' ' show a constant oscillatory character without any presence of dissipative 
effects (see figure I6.9t . The properties of the Hamiltonian constraint violation and L2-norms 
illustrated for the case of an evolution dominated by the F-mode, remains valid also when the 
other overtones are excited and the simulations preserve their stability and absence of dissipation 
also for longer evolutions. 

The simulations carried out with the excitation of the first overtone (HI) show a behaviour 
near the surface that deserves some attention. In figure 16.101 we plot the time profile of the 
two fluid radial perturbations 7( 1 >°) and H^ 1 ' '. The two upper panels have been obtained by 
performing the average (I6.66t for the interior spacetime, while on the bottom the last three grid 
points that are near the surface have been neglected from the average operation. We can first 
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Figure 6. 10: Time evolution of the spatial average do"66l for the radial variables 7^°) and H^'°\ The 
pulsations have been excited with the eigenfunction of -y^ 1 ' ) associated with the first overtone. On the top, 
the averages have been calculated in the whole interior space time. On the bottom, the three grid points 
near the surface have been neglected. 



notice that the enthalpy H^ 1 ' ' preserves the same oscillating properties while the velocity t^ 1 ' ) 
has a higher amplitude and a smoother oscillating dynamics when the points near the surface 
are neglected. By investigating the movie of these simulations for the variable -y^ 1 ' ), we have 
noticed the presence of small spurious oscillations near the surface. This numerical noise is not 
continuous but has a random character, which is sufficient to modify the evolution of this variable. 
In order to have a better understanding we have also analyzed the simulations for the second (H2) 
and third overtones (H3). Similarly to the case of the F-mode, the spatial average profiles of these 
time evolutions do not present any spurious oscillations near the surface. For the second overtone, 
this behaviour is shown in figure 16.1 II However, when we study the corresponding movies we 
notice also in these cases the presence of random and very small spurious oscillations near the 
surface that decrease by increasing the resolution of the meshes. From the analysis of the movies 
and the eigenfunction profiles, we can argue that for the H2 and H3 normal modes the presence of 
a node near the surface seems to prevent the propagation of these oscillations along the star and 
actually reduce their effects. On the other hand for the F-mode, the absence of micro-oscillations 
in the time profiles of figure 16/71 seems more due to the smallness of these numerical oscillations 
with respect to the average value of the physical pulsations, about two order of magnitude less. 
This motivation is confirmed by the analysis of the non-linear simulations when the F-mode will 
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Figure 6.11: Time evolution of the spatial average profiles of the radial variables --y^ 1 - ) and H^- 1 ' '. The 
radial pulsations have been excited with the eigenfunction of ^f 1 °) associated with the second overtone. 



display this noise. As a result, the presence of these oscillations has to be taken into account 
during the implementation of the matching conditions of the non-linear perturbations. 

The spectral properties of the oscillations are studied by performing an FFT of the time pro- 
files. In figure 16.121 we show the power spectrum of seven simulations where every time only 
one of the first seven normal modes has been excited by the initial configuration. Furthermore, 
we have performed a 12 ms simulation for radial pulsations excited by an initial Gaussian pulse. 
The relative spectrum is shown in the lower panel of figure 16.121 where the frequencies of the 
normal modes are labelled with a circle. The radial modes are determined in time domain with an 
accuracy to better than 0.3%, see table 16.21 

Eventually, we determine the convergence rate for the radial evolution which is given in ta- 
ble |0]for the four radial perturbations. 

The amount of pulsation energy contained in the radial perturbations can be determined with 
the expressions of the relativistic kinetic and potential energy derived with a variational analysis 
E2l 12311761 . and lITTll . In particular, according to the initial conditions that we have set up for 
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Figure 6.12: Power spectrum of the radial perturbation r y < ' 1 '°\ which has been determined by an FFTof the 
time profiles. In the upper panel we show in the same plot the spectrum of eight different time evolutions, 
where everytime a single radial mode has been excited. In the lower panel,we show the spectrum of a time 
evolution excited by an initial Gaussian pulse. The excitation of the radial modes is evident in both cases. 
The first eight frequencies are compared to the values determined with a code in the frequency domain and 
are shown with a circle. 



the radial perturbations we can derive the pulsation energy introduced on the initial Cauchy sur- 
face. For an initial vanishing Lagrangian displacement and a specific eigenfunction for the radial 
velocity -yC-' ), the initial pulsation energy is given by the kinetic energy fPTl : 



E^ = 2irW (AW 



Us 



(6.67) 



dr (rV A 7 M) , 

where we have used the relation ( I3.49I ). and jn'°^ is an eigenfunction of -y^ 1 ' ) and A^ 1 ' ) its 
amplitude. In table !6.4l the oscillation energy for the first seven normal modes has been calculated 
with the expression (I6.67t for an amplitude A^ 1 ' ^ = 0.001. 
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Convergence 


7 (1,U) H (l,0) £(1,0) ^(1,0) 




2.04 2.06 2.07 1.53 



Table 6.3: Convergence test for radial perturbations: cr^ 1 ' ) denotes the convergence rate in the L 2 
norm. 



6.4 Linear, axial non-radial oscillations 

The dynamics of axial perturbations on a static star is described by the system of two differential 
equations: the odd-parity master equation (I3.75t for the axial master variable ^C 3 ' 1 ) and the con- 
servation equation ( 13.761 . which is satisfied by the redefined axial velocity perturbation Z^ ' 1 ). In 
section (I3.5t we have already shown that the functional form of the conservation equation (13.76b 
allows us to discern the dynamical degree of freedom of the axial gravitational spacetime from 
the stationary frame dragging profile produced by the presence of an axial differential rotation. 
Hence, the solution of the master equation ( 13.751 can be decomposed in two parts d3.78t . 

= *£2 + ^ 0,1) , (6-68) 

namely the solution VP^2 °f tne homogeneous equation ( 13.791 plus a particular static solution 

f (o,i) 



14, j of equation d3~S0l 
6.4.1 Numerical algorithm 

The homogeneous equation ( 13.791 ) is integrated by means of a standard leapfrog method, which 
is an explicit second order and three level method [91]. By using a centred finite differential 
approximation in space and time we can write the discrete approximation of equation ( 13.791 ) as 
follows: 



-- — A ; 2 - +"j A ; 2 ' • • ■ 2A / ' • "• < 6W > 

Thus, the value of V^om * s updated in time with the following expression: 



A/ 2 ib n — ib n 

= ^1 - ^r 1 + f^ 2 w + i - + ^-1) + °y At 2 yj+i 2A ; j - 1 + Vi At 2 ^ . 

(6.70) 

The coefficient v is the propagation speed of the wave and a, V are coefficients depending only 
on the coordinate r: 

« = e*-\ (6.71) 



(2M 



a = e 29 -Air(p-p)r) , (6.72) 
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y = «»(£-*r»-fl-^]. 



The particular solution is obtained from equation (I3.80I) . which is second-order discretized 
in space and written as a tridiagonal linear system, which is then solved using a standard LU 
decomposition [91 ]. The components of the LU decomposition are given by 

a^7_a + bj1% + c^+i = fj (6.74) 
where the coefficients &j , bj , and dj are 

1 /2m , , \ e 2A 

4vr(p-p)r — - , (6.75) 



Ar 2 V r 2 y 2Ar ' 

b = Ve-™-^T 2 , (6.76) 

1 /2m \ e 2A 

£ = A72 + (jF- 47r (p-p)rJ^, (6.77) 

/ = 167T (4vrpr 2 + ^) e 3 ^ ' 1 ) + leTrre^' 1 ) . (6.78) 
The dragging of the inertial frame can be determined by using expression ( 13.821 ) and the relation 



( I2.96I ). which connects the particular solution ^j, ' 1 ^ with the metric perturbation components 
k 1 ™ and /iq 11 . Alternatively, the metric variable k^ as well as the related frame dragging uj lm can 
also be determined directly by the ordinary differential equation d3.83t . We can apply the same 
procedure used for solving the particular solution ^p 0,1 ^. In this case the LU decomposition 
reads: 

Qj^o, j—i + bjk^j + c/fco, j'+l = fj (6-79) 
where the coefficients dj , bj , and dj are 

e -2A 27rr 

Ar z Ar 
~ , . 1(1 + 1) 4m 2e~ 2A 

b = 87r(p + p)-^ r ^ + ^--— 1 -, (6.81) 
e -2A 27rr 

5 = a^-aF^+p)' (6 - 82) 

/ = levre*^ - 1 ) . (6.83) 



6.4.2 Boundary conditions and initial configuration 

The boundary conditions at the origin, stellar surface and at infinity are implemented in accor- 
dance with the discussion given in section 13.51 At the origin we impose the condition d3.86t at 
the first grid point n = Ar as follows: 



i-t-i 



(6.84) 
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At infinity we use a standard outgoing wave condition. We carry out the integration of the wave 
equation ( 13.791 ) for the homogeneous solution and the ordinary equation (I3.80t for the particular 
solutions on the whole numerical grid without imposing any condition at the stellar surface. We 
have found that the numerical solutions satisfy the continuity of the function t// ' 1 ) and its first 
spatial derivatives t// 0,1 ), which are the requirements prescribed by the junction conditions. For 
the two point boundary value problems (13.801) and ( I3.83I) . we have imposed at the centre and at 
infinity the conditions discussed in section 1331 

The system of equations (I3.75t and d3.76t requires the specification of three functions on the 
initial Cauchy surface, namely 

j(o,i) = ^(0,1)^(0,1)^(0,1)) . (6 . 85) 

In the initial data we specify the original axial perturbations /3( 0,1 ) and we determine the function 
0(0,1) accor di n g to equation (13.77b . 

We can set up two independent classes of initial conditions for the first-order fields: 

1) a stationary differentially rotating star, where the amount of differential rotation can be 
determined by specifying the profile of the axial velocity /?( 0,1 ). The functional form of 
this perturbation can be derived as we see later from one of the theoretical rotation laws 
used in the literature. The profile of the master function ^C 3 ' 1 ) is then obtained by numer- 
ically solving equation d3.80t . Therefore, the set of initial data is given by the following 
expression: 

j(°' 1 ) = (0,0,/3(°' 1 )) , (6.86) 

2) In the second case, we may consider a non-rotating star, Z^ ' 1 ) = 0, and arbitrarily fix the 
metric perturbation ^hom an d its time derivative l I / hom,t> 

4°'M<'i<>) ■ < 6 - 87) 

This initial condition can be used to understand how the second-order metric perturbations 
can be affected by the radial pulsations of the star. In fact, it is well known that the odd- 
parity grativational wave signal arising from a non-rotating compact star can only contain 
the imprint of the spacetime w-modes. The question is then what this particular signal 
looks at second order for a star that is also pulsating radially. In practice, we generate w- 
modes ringing at first order and look at the corresponding signal at second order and at its 
dependence on the radial pulsation of the star. The ui-mode ringing is induced in a standard 
way, i.e. by means of a Gaussian pulse of GWs impinging on the star. 

Let us describe in more detail the choice of the initial conditions for a differentially rotating 
star ( I6.86I ). The differential rotation law of a neutron star is unknown. An accurate description 
of a differentially rotating configuration for newly born neutron stars should come out from the 
numerical simulations of core collapse, references I117ll3*3*ll3*31l . However, a set of rotation laws 
has been introduced in Newtonian analysis, whose main motivations are: mathematical simplicity 
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and the satisfaction of Rayleigh's stability criterion for rotating inviscid fluids: d {w 2 VL) jdw > 
0, where w = rsin# is the cylindrical radial coordinate. The "j-constant law", one of these 
Newtonian laws, has also been extended to the general relativistic approach where the dragging 
of the inertial frame must be taken into account 16411631. 



In this thesis, we provide the initial profile for the odd-parity velocity perturbation pi^. by 
using an expansion in vector harmonics of the velocity perturbation of a slowly differentially 
rotating star. In slow rotation approximation the covariant velocity perturbation is the following: 

SuiQ' 1 ) = (o, Su^A = (0, 0, 0, r 2 sin 2 9 (O - u)) , (6.88) 

where Q = f2(r, 9) is the angular velocity measured by an observer at infinity which describes 
the stellar differential rotation, while the function lo = u;(r,9) denotes the dragging of inertial 
frame associated with the stellar rotation. In barotropic rotating stars, the integrability condition 
of the hydrostatic equilibrium equation requires that the specific angular momentum measured by 
the proper time of the matter is a function of $7 only [63], i.e., u l u^ = In the slow rotation 

case this condition leads to the following expression: 

6j(°'V (Q) = uHuf^ = e- 2 V sin 2 9 (O - u) , (6.89) 

which is valid up to first order in f2. The choice of the functional form of j (O), and hence 
ijj^ ' 1 ) (f2), must satisfy the Rayleigh's stability criterion against axisymmetric disturbances for 
inviscid fluids: 

df/dQ < , (6.90) 
where j~is the specific angular momentum 

f=(p + p)u<f ) /po, (6.91) 

and po the rest mass density. The specific angular momentum f is locally conserved during an 
axisymmetric collapse of perfect fluids 11041 . A common choice for j(Q) that satisfies these 
conditions 1531 loH is the following: 

5j {0 ^(n) = A 2 (n c -n) , (6.92) 

where Q c is the angular velocity at the rotation axis Q(r = 0), and A is a constant parameter 
that governs the amount of differential rotation. Equations ( 16.891 ) and d6.92l > give the following 
expression for the rotation law: 



A 2 Q c + e- 2 *r 2 sin 2 9 u(r, 9) 

e -2$ r 2 s i n 2 



This equation in the Newtonian limit reduces to the j-constant rotation law used in Newtonian 
analysis [49]: 

_ . A 2 n, 
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Figure 6. 13: Wave form in semi-logarithmic and logarithmic scale of the quadrupolar component (I = 2) 
of the axial master function vj/t * 1 ) scaled by the stellar mass M. The excitation of the first w-mode and its 
strongly damped ringing phase are evident for 1.78 < log [t/ (2M)] < 1.95. The late time power-law tail 
is also in accordance to the theoretical results. 

A uniformly rotating configuration with $7 = Q c is attained for high values of A, namely for 
A — > oo. On the other hand for small values of A, the law d6.93t describes in the Newtonian limit 
a configuration with constant angular momentum. The only non-null vector harmonic components 
of 5u a are given by the following expansion: 

i in 

where /3( 01 ) is the scalar function of the coordinate r that has been defined in Eq. d2.70t and that 
can be derived by the inner product with the basis of axial vector harmonics: 

= (Sn^,Sjr) = J BmOdOdtSu^Spi* , (6.96) 

where ~f ab is the contravariant unit metric tensor of the sphere S 2 . 
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Figure 6.14: Wave form in semi-logarithmic and logarithmic scale of the component I = 3 of the axial 
master function vIa ' 1 ) scaled by the stellar mass M. As for the quadrupolar case, the to-mode excitation, 
the ringing phase and the long term time decay is clearly present. 



In order to determine the initial profile of the axial velocity we can introduce equation ( 16.931 1 
into the velocity perturbation (16-881 )- and obtain: 



Su 



(0,1) 



Ah 



A 2 + 



-2$ r 2 gm 2 , 



! sin 2 en c + Y^ k l m S] 



Im 



(6.97) 



lm 



where we have used the relation ( 13.821 that connects the frame dragging function with the metric 
perturbations. It is worth noticing that the axial velocity d6.97t contains as first term the Newto- 
nian j-rotation law up to an exponential factor (hereafter for simplicity, we will call this term the 
nearly Newtonian j-rotation law), while the second part accounts for the frame dragging. When 



we introduce equation d6.97l > into the inner product d6.96b . we can easily determine the expres- 
sion for the nearly Newtonian term while the relativistic correction requires more attention. In 
fact the metric variable k 1 ™ is itself the unknown of the differential equation ( 13.831 and the inner 
product ( 16.961 ) of this term produces quantities which are products of different harmonic indices 
(l,m). In order to decouple the various terms of the relativistic corrections we will assume that 
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Figure 6.15: Power spectrum of the gravitational signal producedby the scattering of an axial gravitational 
wave on a spherical non-rotating star. The solid line refers to the quadrupolar term I = 2 while the dashed 
line to the I = 3 case. The frequencies of the first u>-mode are shown with a triangle for I = 2 and a circle 
for I = 3. The curves are normalized with their maximum values assumed at the peak. 



the dominant contributions will be provided by the k 1 ™ which has the same harmonic index as the 
nearly Newtonian rotation law. When we have determined the harmonic expansion of the nearly 
Newtonian part, we have found a well posed solution only for A > e~®( R ^R s . This relation and 
the values assumed by the metric field <E> in the stellar model considered in this thesis allow us to 
give the following estimation: 



< 



< 1, 



2 - A 2 + e- 2 *r 2 sin 2 o 
which can be introduced in the relativistic term of equation d6.97t to give: 



(6.98) 



Su 



(0,1) 



A 2 r 2 sin 2 6 e-® 
A 2 + e" 2 *r 2 sin 2 9 



l in 



(6.99) 



In i 



Su 



(0,1) 

ao 



where ao G (0.5,1]. The expansion of this law in odd-parity vector harmonic provides the 
following result: 

for I even , 

(6.100) 

I /4 0,1 H° for I odd, 

where the components with odd I are given by the following expression: 



,o(0,l) 

no 



e 9 n c f m (x,A) + a e-' s 'k l , 



(6.101) 
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Figure 6.16: Profiles of the I = 3 component of the axial velocity perturbation ^q' 1 ^ in km, determined 
from a j-constant rotation law with A = 15 km and a T = 10 ms rotation period at the axis. The axial 
velocity associated with the nearly Newtonian j -law is shown with a solid line, while the dotted and dashed 
lines denote the velocity (16. 101 i with a = 0.5 and a = 1 respectively. 



where x = re * (not to be confused with the fluid tortoise coordinate). The functions fa for 
/ = 1 and I = 3 are given by the following expressions: 



f (o,i) 

J10 



f (o,i) 

J 30 



-3.069A 2 
-0.781.4 2 



A 2 



In 



Vx + Va 2 + x 2 



a 2 

1 + 7.5^ 



6A 2 



xV A 2 + x 2 



5A 2 \^ \fx^fWX 



1 + In 



AX 2 J y/^ A 2 + x 2 



X 



(6.102) 
(6.103) 



which have been derived by imposing the condition A > R s e~®( Rs \ 

In the limit of A — > oo, the function fa (r, A) behaves correctly, i.e. the nearly Newtonian 
part assumes the following expressions: 



f e~*0,r 2 sin 2 9 for I = 1 



lim e n c / zo 5, 

A— >oo 



10 



(6.104) 







for I > 3 . 



which corresponds to a uniform rotating configuration, where Q, = Q c is the angular velocity 
measured at infinity. The aim in this work is the analysis of the gravitational signal produced by 
the non-linear coupling. Therefore, we will consider only the case I = 3 as the dipolar term I = 1 
does not produce gravitational waves. 

In figure l6~T6l is plotted the axial velocity perturbation ft^ 1 ^ . The figure shows that in the 
two extreme values of the constant an, the solutions disagree up to ten percent. In the following 
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Figure 6.17: Profiles of the I 



3 components of the axial fluid perturbation 0^ m ^ m 1 



(solid 



line) and its spatial derivative P^' 1 ^ in km 2 (dashed lines), determined for a j -constant rotation law with 
A = 15 km and a T = 10 ms rotation period at the axis, and with a = 0. 



simulations we will use for simplicity only the nearly Newtonian term and we obtain results which 
are correct to better than ten percent. 



6.4.3 Simulations for axial non-radial perturbations 

In this section, we carried out numerical simulations of linear axial perturbations for the two 
initial configurations described in section l6".4.2l 

There are no fluid oscillations associated with linear axial perturbations, therefore these are 
characterized by pure spacetime oscillations represented in terms of a set of high frequency and 
strongly damped w-modes. These spacetime modes can be excited dynamically by studying the 
scattering of an axial gravitational wave on a compact star. This is a standard procedure exten- 
sively used in the literature llT2ll36l . The gravitational wave can be modelled by an impinging 
initial Gaussian profile, 



(0,1) 
horn 



A (0,l) e -q(r-r ) 2 



t=t 



(6.105) 



where A^ ' 1 ) designs the initial wave amplitude, while the constants q and tq control the wide 
and initial position of the Gaussian respectively. We set up two simulations for the / = 2 
and / = 3 components of the axial perturbative fields. The initial Gaussian pulse is at 40 km 
far from the centre of the star, has an amplitude A^ ' 1 ^ = 0.1 km and q = 1.25 km~ 2 . The 
wave forms determined by an observer at 150 km are shown in figures (I6.13t and (16.141 . The 
scattered signal in both cases exhibits the characteristic excitation of the first w-mode at about 
log[t/(2M)] = 1.8 evolution time, followed by the ringing phase which is strongly damped by 
the emission of gravitational radiation. The ringing of the QNM is more evident in the lower 
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Figure 6.18: The upper panel displays the stationary axial master function ^p ' 1 "*, in km, for a nearly 
Newtonian j -constant rotation law with A — 15 km and a period T = 10 ms at the rotation axis. The 
solution of equation J3.80l > is shown as a solid line while the dashed line is the solution found indirectly by 
first solving equation d3.83i for the variable fcg°' 1 ' ) and then using the definition J6. 108i . 



panels of figures (16.131) and (16.141) where the signal is plotted on a logarithmic scale. After this 
phase, the master function ^h'^ shows the typical decreasing behaviour predicted by the theo- 
retical power law Vpj^ ~ £~( 2 '+ 3 ) ||26ll . For an evolution time of 8 ms a linear regression of the 
tail of the signal provides the following values: 

(o i) i t~ 7m for I = 2 . 

= fOTl = 3 . < 6 '°« 

The power law gets closer to the theoretical value for longer evolutions. 

The excitation of the first w -modes for the two harmonic indices I = 2, 3 is confirmed by the 
analysis of the spectral properties of the signal. We perform an FFT of the signal part that starts 
with the excitation of the spacetime mode. The results in figure l6~T31 show the characteristic broad 
shape of w-modes which is due to the short values of their damping time, respectively r = 29.5 fj,s 
and r = 25.2 us for the I = 2 and I = 3 cases. The curve peaks at the v = 10.501 kHz for the 
quadrupolar case and at v = 16.092 kHz for the I = 3 case, which are denoted in figure l6~T51 wifh 
a triangle and circle respectively. For these two harmonic indices these values are the frequencies 
of the first curvature w-mode, which have been determined by Gualtieri's numerical code in a 
frequency domain approach l47l . 
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Figure 6. 19: The upper panel displays the I = 3 component of the axial velocity perturbation ft ' 1 ', in 
km, for a nearly Newtonian j-constant rotation law with a period T = 10 ms at the rotation axis and for 
different values of the differential parameter A (in km). The associated solutions of equation J3.80I of the 
stationary master function tpp ' 1 ^, in km, are shown in the lower panel. 



The other initial configuration for the linear axial perturbation describes a stationary differen- 
tial rotation induced by the axial velocity perturbation. In particular, the structure of this pertur- 
bative framework allows us to investigate the (l,m) component of the axial perturbation ifi^ ' 1 ', 
which is induced by the (/, m) component of the axial velocity perturbation /?( 0,1 ). As illustrated 
in section lrj.4.21 we are going to implement an axial velocity perturbation described by the func- 
tion d6.101b . which was derived from the relativistic j-constant rotation law. The first point to 
clarify in equation d6.101t is the amount of the relativistic correction due to its second term with 
respect to the nearly Newtonian j-constant law. This issue can be studied with equation d3.83t for 
the metric gauge invariant quantity k^' 1 ^. We have introduced the functional dependence (I6.101t 
for Z^ ' 1 ) into the source term of equation (I3.83t and solved the equation first for a nearly Newto- 
nian rotation law, i.e. = 0, and secondly for «o = 0.5 and ao = 1, which are the two extreme 
values of the relativistic correction to the I = 3 velocity component. In addition, we choose the 
following value of the differential rotation parameter A = 15 km and a T = 10 ms rotation 
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period at the rotation axis of the star. The choice of A is motivated by the regularity of the I = 3 
component of the axial velocity (see section l6.4.2t . that leads to the condition A > e®( Rs ^ R s . 
For higher values of A, the harmonic components (5^^ with / > 1 will decrease, as the rotation 
law approaches the uniform rotation configuration. On the other hand, the angular velocity Q c is 
the physical perturbative parameter that controls the strength of the axial perturbations. Its values 
must satisfy the requirement imposed by the "slow rotation approximation", i.e. the dimension- 
less perturbative parameter e = Q c /Qk *C 1, where 0,% is the Keplerian angular rotation which 
describes the mass-sheeding limit of the stellar model under consideration. In order to have a 
simple extimate, we consider the Qk of a uniformly rotating star, which can be extimated with 
the empirical formula [50 38 1 which is accurate to better than 10 percent. This formula is based 
on the classical value up to a corrrective factor and is given by the following expression: 



n K = 0.625 y/M/R* , (6. 107) 

where M and R s are the mass and radius of the non-rotating star in hydrostatic equilibrium. 
For our stellar model we get Qk = 0.0324 km' 1 which corresponds to a rotational period 
Tk = 193.996 km = 0.64 ms. Therefore, the dimensionless parameter is e = 6.45 x 10~ 2 . 

The numerical integrations of equation (13.831) with the method explained in section l6".4.1l are 
shown in figure l6~T6l The two solutions obtained for «o = 0.5 and qo = 1 disagree up to 5 
percent, while the Newtonian differential rotation law is accurate with respect to the relativistic 
ao = 1 case to better than 10 percent. We will use the nearly Newtonian rotation law in this 
and the next sections, being aware that the results could be accurate within ten percent due to the 
relativistic corrections of the dragging of inertial frame. 

The particular solution of the axial master function ^p ' 1 ^ can then be determined by equa- 
tion (I3.80t . which is shown in figure l6"T%l for the same parameter of the j-constant rotation law 



used above. This solution can also be determined indirectly by first solving equation (I3.83t for 
the metric variable fco an d then getting ^p ' 1 ^ through the definition (I3.74t for the stationary case: 

*(?.!) = ( 2 Jfe< 0>1) - rk ( ff) e-(* +A ) . (6.108) 

The indirect solution, which is shown in figure 16.181 as a dashed line, reproduces the solution 
obtained with the direct method with a maximum error less than 2.3 percent. In the lower panel 
of figure |6~T%1 we show the I = 3 component of the frame dragging function W30, which has been 



found from equation (I3.82t . 

So far, we have used a rotation law with a rotation period at the axis fixed at T = 10 ms. 
However, the linearity of equation (I3.80t allows us to determine the solutions ^p ' 1 ^ with a simple 
rescaling. Let Vl/l ' be a solution related to a differential rotation period T\. The solution 

Ti 

corresponding to a rotational period T2 is given by: 



^(0,1) = ^ ^(0,1) 



V 



T 2 Ti 



2 



V 



(6.109) 
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In figure 16.191 we show the effects of the differential parameter A on the profile of the axial 
velocity /?( 0,1 ) for ao = (I6.101t and the particular solution of the master function ^'^ > ' 1 \ de- 
termined by equation (I3.80t . The I = 3 axial velocity and master function decrease for higher 
values of A, when the star tends to a uniform rotational configuration and the only non- vanishing 
component is the / = 1. In addition, we have noticed that for A > 100 km the axial velocity is 
not accurately described by the expressions (I6.102t and (16. 103b . as it has an irregular behaviour 
near the origin due to appearence of high peaks. 

The rotation energy associated with the differential rotation can be determined with the fol- 
lowing equation [53]: 

£(o,i) * f ° dr [ d9 2Trr A sm6 3 (p + p)e A -*n(n-uj) . (6.110) 
2 Jo Jo 

By neglecting the dragging of the inertial frame u we can determine an upper limit of the rota- 
tional energy [ 53 ] : 

E^ ] = \ [ S dr [ d9 2vrr 4 sin 9 3 (p + p) e^Q 2 , (6.111) 
2 Jo Jo 



where E^ ' 1 ) < E^\ Furthermore, this expression can be expanded in tensor harmonic as 



J rot 

follows: 



i> 



where we have applied the following relation f^ m = r~ 2 e*/3j° , which is valid when the frame 
dragging function uj is neglected. The quantities Qj m are the harmonic components of the angular 
velocity, 

n(r,0) = £n, m Sf m . (6.113) 



J>1 



The upper limit of the Z = 1, 3 component of the rotational energy for a differential rotation with 

A = 15 km and T = 10 ms is e£J$ = 7.90 x 10~ 5 km and = 8 - 62 x km 

respectively. 

The solutions and ko exhibit a second order convergence, while ^p ' 1 '* manifests a 

convergence of first order. This lower convergence rate is due to the discontinuity of VPp^r at tne 
stellar surface. 



6.5 Non-linear axial oscillations 

The dynamical properties of the coupling between the radial and axial non-radial oscillations 
are described by the solutions of the two inhomogeneous partial differential equations (14.62b 
and ( 14.63b . where the source terms are present only in the interior spacetime. The initial values 
adopted for the linear perturbations are able to describe the non-linear coupling for the following 
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configurations: i) a radially pulsating and differentially rotating star, ii) scattering of a gravita- 
tional wave on a radially oscillating star. In particular in the former case, we will see that the 
coupling between the stationary axial velocity with the radial pulsations produces an oscillating 
Ae corrective term of the redefined axial velocity fit 1 ' 1 '. 

6.5.1 Numerical algorithms 

The conservation equation (14.631) for the perturbation ft 1 ' 1 ' exists only inside the star and is 
integrated with an up-wind algorithm. The numerical discretization is given by the following 
expression: 

= 4" + At e*' (6. 1 14) 

where (S^)^ is the second order discrete approximation of the source term given in appendix (fPl. 
Thus, the first and second order derivatives that are present in (S^) . are approximated by second 
order centered finite difference approximations in the internal grid points and by second order 
one-sided finite approximations at the origin and stellar surface. 

The integration domain of the axial master equation (I4.62t is the entire spacetime, where the 
source £^ is present only in the interior spacetime. After several tests, we have found that the 
simulations of the master function vjn 1 ' 1 ) are more accurate when we implement two different 
numerical methods. An Up-wind algorithm for the interior and a Leapfrog for the exterior, for 
studying the coupling between the radial pulsations and differential rotation, and a Leapfrog on 
the whole spacetime for the scattering of an axial gravitational wave on a radially oscillating star. 
The implementation of these two methods is due to the different properties of the source terms 
and the junction conditions (I4.70t for the two cases mentioned above. 

In case of coupling between the radial pulsations and axial differential rotation, we prefer 
to separate the numerical integration of the axial master equation in the interior and exterior 
spacetime. Furthermore, in order to reduce the numerical noise caused by the discontinuity of 
the sources at the stellar surface, we implement a first order accurate numerical scheme in the 
interior spacetime. Therefore, the axial master equation is simulated with two different numerical 
schemes inside and ouside the star. In the interior, the second order PDE H4.62I is transformed 
in a system of first order PDEs which will be integrated with a generalization of the Up-wind 
method [67|. In the exterior, the Regge-Wheeler equation at order Ae is instead updated in time 
by using a second order Leapfrog method, where we extract the values of vJ/C 1 ' 1 ) and SSf;r on the 
surface with the matching conditions. Let us first describe the interior spacetime. We can define 
two new quantities as follows: 

W = ^ini 1) ' Ui=W t , u 2 =W jr ., (6.115) 

where the variable w has been introduce to reduce the number of indices in the discrete equations. 
With the definition of the vector u = [ui, 112 ] T we can write the master wave equation (I4.62t in 
the interior spacetime in the following conservative form: 



u, t + F, r =S 



(6.116) 
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where the flux is given by 

and A is the following 2x2 matrix: 



F = Au, 



(6.117) 



v z gw 

1 



The quantity v gw = e is the propagation velocity of the gravitational signal, and the source S 
is a two dimensional vector with the following components: 



Si 



2m \ 



-e 2 * ( 4vr (p - p) r + u 2 + Vw + 16tt (4vrpr 2 - ™) e 2 *+A£(l,l) 



+16vrre 2,I> - A /3^ 1 )+S v ,, 
S 2 = 0, 



(6.118) 
(6.119) 



where £^ are the sources expressed in equation (ID. II) and are second order accurate in space. The 
variable u is updated at every time step by the following differentiating scheme, 



u 



ra+l 



£ A+ (u« - u^O - ^ A" (u? +1 - u«) + At S? , (6.120) 



where the matrices A + and A are given by, 



Ax 



-v 

V 



V V 
1 V 



The value of the ip^P is then obtained from the definition ( 16.1151 ): 



w 



n+l 



n At 

w H 

^ 2 



+(111)! 



(6.121) 



In the exterior spacetime the system of equation reduces to the Regge-Wheeler equation for the 
master function xJa 1 ' 1 ). This axial wave-like equation is solved with a Leapfrog algorithm as 
well as we did for the first order axial master equation. We do not write the explicit discrete 
expression for the variable vl^ 1 ' 1 ) here as it can be determined directly by equation ( I6.70I ). where 
the coefficients are now given by 



v 
a 

V 



J2$ 



2M 



2* 



6M 1(1 + 1) 



(6.122) 
(6.123) 

(6.124) 



The scattering of a gravitational wave on a radially pulsating star is instead studied with 
a Lepfrog algorithm in the whole spacetime. When we discuss the numerical implementation of 
the junction conditions for non-linear perturbations (section lo.5.2t . we will see that the movement 
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of the stellar surface in a Eulerian gauge leads us to adopt some approximate treatments of the 
matching conditions. As already explained in section 14.2.21 we have chosen in this work to 
impose the junction conditions on the surface T>j c which is always well inside the star, even 
during the phases of maximal contraction. Consequently, the axial junction conditions reduce 
to the continuity of vjK 1 - 1 ) and its derivatives ( 14. 701 ) . We can then integrate equation d4.62t by 
using the explicit discrete scheme (I6.70t . where we introduce in addition the second order finite 
approximation of the source in the interior spacetime. We have noticed that the junction 
conditions are automaticaly satisfied by the simulations. 

6.5.2 Boundary and initial conditions 

The three non-linear coupling perturbations vl/C 1 ' 1 ), v]/^ 1 ' 1 ) and /J^ 1 ' 1 ) are not all independent on 
the initial Cauchy surface, because of the presence of not vanishing source terms in the non-linear 
perturbative equations. A correct initial value for the axial master function 'J/t 1 ' 1 ) must come 
out from the solution of the coupling axial constraints, as in case of first order axial and polar 
perturbations in presence of sources ll82ll . However, since we are more interested in the part of 
the gravitational wave signal which is driven by the stationary radially oscillating sources, we set 
a vanishing value for all three coupling perturbations at t = without solving the constraints. 
As a result, we will see later that after an initial burst of u>-mode gravitational waves due to the 
violation of the constraints the solution relax to the correct periodic solution. 

Let us now discuss the boundary conditions. At the origin we implement the conditions (14.68b 
by applying the same procedure described for the linear perturbations. Hence, we have for the 
metric variable: 

/ \ 2+1 / \ l+l / \ I 

Wi=W 2 (^J , (*l)l = (*l) 2 (j^J , (U 2 )i = (u 2 ) 2 f , 

(6.125) 

and for the fluid quantity 

^) = ^)^y +i . (6 . 126) 

At the outermost grid point we implement the outgoing Sommerfeld boundary condition, where 
the outermost grid point is set far enough from the star in order to not corrupt the gravitational sig- 
nal extraction with the numerical noise coming from the numerical reflection at the grid buondary. 

The internal and external spacetimes are separated by the stellar surface, where the matching 
conditions have to be imposed in order to connect the physical and geometrical descriptions be- 
tween the interior and exterior spacetime. In an Eulerian gauge, the perturbed surface could not 
coincide with the surface of the equilibrium configuration. Therefore, when we compare near the 
surface the perturbations and background tensor fields related to a same physical quantities some 
problems can arise. Let us for instance consider radial pulsations of a static star, and a point xp, 
which is in the region between the static and perturbed surface during a contraction phase of the 
star x(t) < xp < R s , where x(t) = R s + A^ 1 ' ) and is the Lagrangian displacement of 
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radial perturbations. During this dynamical phase xp is outside the star, then the mass energy 
density p should vanish there. However, when we determine at xp the Eulerian perturbation 



where p is the background density, we find this unphysical value p^- ) = p < 0. At any perturba- 
tive order, the dynamical properties of the perturbation fields depend on the quantities which have 
been determined at previous perturbative orders. Therefore, the introduction of this negative val- 
ues for the density into the Ae perturbative equations could leads to an unphysical description of 
the non-linear perturbations. This issue can be solved with different methods [ 101 ]. For instance, 
one may add to the stellar model an atmosphere of low density, which extends beyond the size of 
the neutron star. Otherwise, one can adopt a Lagrangian gauge for the perturbative description. 
In this case, the perturbed and background surfaces concide and the energy density would get the 
correct value 



where the vanishing value is taken at the background stellar surface. Alternatively, we can impose 
the junction conditions on the hypersurface d4.55t . which remains always inside the star. With this 
method we remove the outer layers of the star, but since the density in this region is extremely 
low the mass neglected is less than one percent. Therefore, the wave forms and spectra of the 
gravitational signal are still well approximated by this procedure. The position of the junction 
hypersurface Sj c is evaluated as follows: i) we evolve the radial perturbations and determine the 
maximum movement of the surface with the Lagrangian displacement £^ . The values for the 
particular initial conditions used in our evolution are written in table 16.41 ii) Then, we place the 
surface £ JC at the first grid point that remains inside the star during the evolution. It is worth 
remarking that for an r-grid of dimension J r = 200, the values of C^/ m table 16.41 leads to 
neglect only one point, i.e. the junction conditions can be imposed at grid point Jj c = 199. 
However, since the source of the axial master equation ( 14.621 ) has its maximal amplitude at the 
stellar surface, this removal of a grid point induces a error of about five percent in the gravitational 
signal. 

For the coupling between radial pulsation and differential rotation, we impose the continuity of 
the metric variable vJ/C 1 ' 1 ) and the relation ( 14.701 ) for ^r' 1 ^. For the latter condition we discretize 
the sided external derivative with a first order approximation in order to get the appropriate value 
of the master function vIa 1 * 1 ) in the first point outside the junction surface Ej c : 



P 



,(i,o) 



= P-P 



(6.127) 



Ap m =p (i,o) + £ p > 



(6.128) 




(6.129) 



(6.130) 



(6.131) 
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F HI H2 H3 H4 H5 H6 

E ( n' 0) [10~ 8 km] 35.9 4.2 1.37 0.62 0.34 0.21 0.14 
^J.' 0) [m] 12.65 4.02 2.66 2.02 1.64 1.38 1.19 

Table 6.4: Energy Eq 1 ' ^ and maximum surface displacement £^ of radial pulsations. These values 
correspond to the initial conditions J6.58i and (I6.65> with A^ 1,0 > = 0.001. Every evolution has been 
excited with a single radial mode. 

6.5.3 Coupling between radial pulsations and axial differential rotation 

In this section, we present the results of the non-linear perturbations which describe the coupling 
between a stationary differential rotation and radial pulsations. For the particular choice of the 
initial rotating configuration the non-linear perturbations studied in this section are relative to the 
harmonic index I = 3. The axial rotation of the fluid is described with a good approximation 
by the nearly Newtonian j -constant rotation law which has been derived in section 16.4.21 As 
in section I6.4.2I and [6.4.31 we have specified the differential parameter (A = 15 km) and the 
angular velocity at the rotation axis (f2 c = 2.09 10 -3 km -1 ). This value, which corresponds 
to a 10 ms rotation period at the axis, is small with respect to the mass shedding limit rotation 
rate Q. c /Q,k = 6.4510~ 2 . The stationary profile of the axial velocity perturbation f3^ 0,1 ' for the 
harmonic index I = 3 is given by the expression (I6.101I) . and has been plotted in the figures (I6.16t 
and d6T9T 

The radial pulsations are instead excited for any evolution at t = by selecting a single normal 
mode given by equation d6.58l >. where the initial radial velocity perturbation has the form of 
equation ( I6.65t with an initial amplitude A( 1,0 > = 0.001. The radial simulations related to these 
initial conditions provide the maximum Lagrangian displacements shown in table !6.4l which have 
been determined at the static surface of the star. From these values, we can argue that the region of 
the spacetime where the perturbed surface of the star moves is really confined in a narrow region 
around the static equilibrium configuration. Therefore, the issues related to the negative values of 
the mass-energy density in a Eulerian description can be approximatively described by matching 
the non-linear perturbations on a surface always contained into the star. For more details see 
section l6".5.2l 

We start by investigating the evolution of the axial master function xJa 1 ' 1 ), when the radial 
perturbations oscillate at the frequency of the first overtone HI. In figure l6~20l we show the wave 
form for an observer at 100 km from the stellar centre. This curve presents a first excitation of a 
typical w-mode followed by a monochromatic periodic signal. This interpretation is confirmed by 
the spectrum associated to this signal. Indeed the left panel of figure loTTl displays the broad peak 
of the w-mode, which is centered around the frequency v = 16.092 kHz, and a narrow peak, 
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Figure 6.20: Wave form of the axial master function ^/A 1 ' 1 ) scaled by the stellar mass M, which describes 
the coupling between the linear radial pulsations and the axial differential rotation of a neutron star. The 
rotation is given by a nearly Newtonian j -constant rotation law with a T = 10 ms period at the rotation 
axis and A = 15 km, while the radial pulsating dynamics has been excited with the overtone HI. The 
signal manifests an initial excitation of the first I — 3 w-mode, which is followed by a periodic oscillation 
driven by the pulsating source terms. 



which has the same frequency of the first overtone HI that has been mirrored in the gravitational 
signal at second order. In the right panel, the FFT of the signal has been performed for a higher 
number of time cycles where the signal is dominated by the periodic evolution of the source 
terms. As a result, the energy contained at the frequency of HI is much higher than that of the 
curvature mode. The presence of a periodic dynamics which reflects the radial pulsations into 
these non-linear coupling arises from the structure of the source terms. In fact for a stationary 
axial perturbation, we can schematically write the sources of equations J4.62t and (I4.63t as 

S = J2ln R (*, r) J* it, r) = J2ln R (r) J* (r) e*"»« (6.132) 

a n 

where the superscripts R and NR denote respectively the linear radial and non-radial perturba- 
tions and Lo a are the discrete set of radial modes. It is then evident that the dynamics is entirely 
sustained by the normal modes of the radial perturbations and that the Ae gravitational signal 
\E'( 1 > 1 ) and the axial velocity perturbation fit 1 ' 1 ' oscillate at the frequencies of the radial modes. 
This behaviour is expected from the symmetries of the stellar model. In fact, the radial frequencies 
could be corrected by rotation starting from the perturbative order Ae 2 , as they must be invariant 
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Figure 6.21: Normalized power spectrum of the wave form shown in figure l6"20l The left panel displays 
the presence of the HI radial mode and the I = 3 w-mode in the gravitational signal. In the right panel the 
FFT has been performed for a longer evolution time, thus the energy contained in the HI peak becomes 
dominant and the spacetime mode is not visible. 



for an inversion of the rotation direction Q — > — Q. 

On the other hand, the presence of the first u>-mode seems to have a different origin. This space- 
time mode is only excited during the early phases of the dynamical evolution. As explained in 
section lo". 5 .21 this gravitational burst is due to the initial transient in tlK 1 ' 1 ), which is produced by 
the initial constraint violation. In order to confirm that a transient can generate a burst of grav- 
itational wave characterized by the excitation of a u>-mode, we modify the source of the axial 
master equation ( I4.62t with an artificial transient described by a delta function placed in a point 
inside the star at a given time. This delta function is then immediatly removed after a time step. 
The wave form and the spectrum of the resulting simulation carried out for the harmonic indices 
I = 2, 3 are shown in figure l6~22l where the excitation of the first u>-mode is evident. 

When we extend this analysis to simulations where the radial pulsations are excited with 
higher overtones and with the fundamental mode, the wave forms and the spectra have the same 
features described for the HI case. However, an interesting amplification is noticed in the grav- 
itational signal when the radial perturbations pulsate at frequencies close to the axial w-mode. 
In figure |6~2"31 we show with the same scale the axial master function for six simulations, 

where the initial data for the radial perturbations are provided each time by one of the first six 
overtones. It is evident that the function $( 1 > 1 ) increases in amplitude when the order of the radial 
normal mode increases from the first to the fourth overtone, while the amplitude decreases for the 
fifth and sixth overtone. The amplitude of vIH 1 ' 1 ) related to the radial H4 evolution results is about 
sixteen times the amplitude of vl/t 1 ' 1 ) for a HI evolution. For this stellar model, the spacetime 
mode for the harmonic index I = 3 has frequency 16.092 kHz which is between the frequencies 
of the third and fourth overtones, 13.545 kHz and 16.706 kHz respectively. In addition, we 
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Figure 6.22: The wave forms (left panel) and spectra (right panel) of the master function tp, which 
corresponds to the toy model studied in section l6.5.3l The presence of an arbitrary transient in the source 
of equation J6.5.3> . which has been induced by a delta function, excites the first w-mode. 



can notice that this effect is present although the energy contained in the radial pulsations and 
the maximum displacement of the surface decreases proportionally for higher radial modes (see 
table 16.41 . Considering the structure of the axial master equation (I4.62L this amplification is 
certainly due to a resonance between the axial potential V, which contains the properties of the 
QNM of the spacetime, and the source that is pulsating at the radial eigenfrequencies. In fact, 
we have noticed that this effect disappears when the axial potential is arbitrarily removed from 
the equation. In addition, the fluid variable /^K 1 ' 1 ) that obeys the conservation equation (14.631 
does not show amplification and decreases proportionally with the order of the radial mode. In 
analogy with a forced oscillator, the amplification of the signal appears when one of the natural 
frequencies, determined by the form of the axial potential V, is sufficiently close to the frequen- 
cies associated with the external force term. Moreover, it is interesting to notice in figure l6"2"3*l that 
in accordance with the broad spectral peak of the w-modes (figure l6~T5l . this resonance affects a 
wide frequency band around the w-mode. 

Consistently with the previous results, the power emitted at infinity by gravitational radiation ex- 
hibits the amplification present in the axial master function. In figure l6~24l we determine the rate 
of the radiate energy for the six simulations which have been described before and were shown 
in figure l6~23l The power displays two distinct contributions: i) a first peak which is due to the 
energy emitted by the first u>-mode which comes out from the initial constraint violation, ii) the 
periodic emission due to the pulsations of the source. However, the relative strengh of these two 
effects changes when we excite the radial oscillations with different eigenmodes. In fact, the pe- 
riodic emission first increases proportionally with the order of the radial overtone excited (HI, 
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Figure 6.23: Comparison of six wave forms of the axial master function tp^ 1 ' 1 ' given in km, where the 
axial perturbations are described by the same differentially rotating configuration illustrated in figure lo*?20l 
On the other hand, the radial pulsations have been excited each time with one of the six radial overtone 
(Hl-6). The function t// 11 ', which is plotted on the same scale, shows a resonance effect when the radial 
perturbations pulsate at H4 overtone, whose frequency is close to the first lomode. 



H2 ,H3). In correspondence with the H4 radial overtone it reaches its maximum and dominate 
completely the radiation, and eventually it decreases for the H5 and H6 overtones. 

Now, we can study the coupling between differential rotation and radial pulsations that os- 
cillate at the fundamental frequency F-mode. From the previous analyses, we can expect a grav- 
itational signal with at initial excitation of the axial u>-mode, which is immediately followed 
by periodic pulsations with the F-mode frequency which are driven by the source. In addition, 
the axial potential V should lower the second part of the signal as the F-mode frequency is far 
from the first u>-mode. This expected behaviour is present in the wave forms shown in the upper 
panels of figure 16.261 Unfortunately, the periodic part of the signal exhibits also some micro- 
oscillations. They are the spurious numerical oscillations discussed in section 16.3.41 which arise 
from the regions at low density near the stellar surface and that seem to depend on the profile of 
the eigenfunction. The signal can be partially improved by neglecting the outer layers of the star. 
In the left column of ngure l6~26l the junction conditions have been imposed on the hypersurface at 
radial coordinate r = 8.64 km, which corresponds to neglect the 0.3 % of the stellar gravitational 
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Figure 6.24: For the same six simulations described in figure I6T231 we plot the power emitted in gravita- 
tional waves at infinity. The six panels do not have the same scale. 



mass. The waveform (upper panel) and the spectrum (lower panel) clearly diplay the corruption 
of the signal due to the numerical oscillations, which appear at the frequencies of the radial mode 
overtones. In the right column, the matching surface is at r = 7.75 km and the mass neglected is 
the 6.3 % of the total mass of the star. Despite the lower amplitude of the oscillations we can now 
notice that the waveform and spectrum are closer to the expected form and the numerical noise is 
reduced. In the spectrum, we can distinguish the excitation of the F-mode and / = 3 w-mode as 
well as some small peaks due to the numerical noise, especially at the frequency v = 13.4 kHz. 
In order to estimate the error introduced by the micro-oscillations in the "expected" wave forms, 
we first proceed to eliminate the first part of the signal characterized by the initial transient. Then 
we interpolate the periodic part of the wave form with a trigonometric function, 

f nt = c sm(2irv F t + c 1 ) (6.133) 

where vp is the F mode radial frequency and cq and c% are two free parameters, which are de- 
termined by the interpolation algorithm. It is worth noticing that the interpolation functions for 
the two signals shown in figure l6~26l will be in general different, since the amplitude of the wave 
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Figure 6.25: Time evolution of the second order fluid perturbation /j^ 1 ' 1 ) in fcm -1 , which is related to the 
axial velocity through the definition (14.651 . The perturbation /j^ 1,1 ) has been averaged at each time step on 
the interior spacetime with the formula J6.66I . The curves plotted in the six panels refer to six evolutions 
where the radial pulsations have been excited every time with a single overtone of the radial modes. In the 
upper panel of the^zrsf column, we compare the arising from the HI (solid line) and H2 (dashed 

line) radial pulsations. With the exception of the top panel of the first column, all the figures have the same 
scale. 



form depends on the position of the matching surface. Thus, we separately estimate the devia- 
tion between the evolved wave form and its interpolation function f mt by using the following 
expression: 



< A^ >= 



1 J 

-Y 



\ 2 

fmt I 
J 3 J ' 



(6.134) 



where J is the number of grid points. Equation (I6.134t provides < A^ >= 1.76 x 10 -8 km and 
< A^ >= 4.45 x 10 -9 km for the wave forms extracted on the junction surface at r = 8.64 km 
and r = 7.75 km respectively. The corruption of the signal is then reduced in average by a factor 
of four at r = 7.75 km. This improvement can be also noticed in figure |6~T71 where the wave 
forms and the interpolation functions are shown for the junction conditions at r = 8.64 km (top 
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Figure 6.26: Wave forms and spectra of the axial master function vIa 1 ' 1 ), in km, for the coupling between 
differential rotation and radial pulsations excited by the F-mode. The left column displays the wave form 
(top) and spectrum (bottom) for a simulation where the junction conditions have been imposed on a hy- 
persurface at r = 8.64 km. The right column shows the same quantities, but now the matching surface 
is at r = 7.75 km. In accordance to the results shown in figure lo"23l the excitation of the F-mode and 
w-mode for I — 3 are evident. In addition, the curves manifest the presence of spurious micro-oscillations 
that appear at the frequencies of the radial higher overtones (vertical dotted lines). This noise is reduced 
when the matching is imposed on a more internal hypersurface. 



panel) and r = 7.75 km (middle panel), while the deviation A^ 1 ' 1 ) = vJ/C-' 1 ) — f int is shown in 
the lower panel. 

The perturbative approach used in this thesis does not include the back reaction analysis, 
i.e. the damping of the radial oscillations or the slowing of the stellar rotation due to the energy 
emitted by the Ae gravitational radiation. This effect could be studied by investigating the third 
perturbative order, which is beyond the aims of this thesis. However, we can provide a rough 
extimate of the damping time of the radial pulsations. We can assume that the energy emitted 
in gravitational waves is completely supplied by the first order radial oscillations, and that the 
power radiated in gravitational waves is constant in time. Thus the damping time is given by the 
following expression: 

p(1.0) 

J 1 ' 1 ) = ^ n (6 135) 

" £(l,i) . ' 

where En i s the energy of the radial pulsations relative to the specific initial eigenf unction (see 
table 16.41 ). while < > is the time averaged value of the power emitted by the Ae axial 
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Figure 6.27: The axial master function ^ in km, (dashed line) and interpolation function f mt (solid 
line) are plotted for the coupling between differential rotation and radial pulsations excited by the F-mode. 
In the upper panel the curves are obtained for a simulation where the junction conditions have been im- 
posed on a hypersurface at r = 8.64 km. In the middle panel, the same quantities are now determined 
by matching at r = 7.75 km. In the lower panel, the difference A^ 1,1 ) = yjK 1 ' 1 ) — f mt is plotted for 
r = 8.64 km (dashed line) and r = 7.75 km (solid line). 



gravitational wave, which has been calculated by averaging equation (14.67b . The calculations 
carried out with the initial data described before for the linear perturbations provide the results 
reported in table 16.51 From the definition (I6.135t and the bilinear dependence of the non-linear 
perturbations Ae, the damping time depends only on the non-radial parameter e, as ~ 
e -2 . In the third row of table !6.51 we write the damping of the radial pulsations due to the coupling 
in terms of oscillation periods for any radial mode, namely 

N osc = , (6.136) 

Ml 

where P n = v~ x and v n is the eigenfrequency of the radial normal mode. 

The Ae non-linear perturbations are bilinear with respect to their perturbative parameters. In 
order to test this property in the numerical simulations, we have performed a first simulation by 
fixing a basic rotation period for the axial perturbations To = 10 ms and an initial amplitude 
of the velocity perturbation Aq = 0.001, and we have determined the averaged values of the 
non-linear perturbations ^q 1 ' 1 ^ and (3q\ Then, we have carried out several simulations with 
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HI 


H2 


H3 


H4 


H5 


H6 



[io- 13 ] 1.85 x 1(T 6 6.83 x 1(T 2 4.85 56.03 157.01 19.08 4.45 
r^' 1} [ms] 6.49 x 10 9 20.49 x 10 3 94.15 3.69 0.72 3.67 10.51 
N osc 1.39 x 10 10 1.408 x 10 5 971.58 49.99 12.07 72.7 8 240.77 

Table 6.5: Rate of energy E^ 1 ' 1 ) emitted in gravitational waves at infinity with the coupling between 
radial pulsations and axial differential rotation. Estimated values of the damping times and number 
of oscillation periods iV osc necessary for the non-linear oscillations to radiate the whole energy initially 
contained by the radial modes. 



different periods T n = uTq, and radial amplitude A m = mA$, where n,m £ N. The associated 
values for VP^ = fe\Pg^ and (5^^ = kfi^' 1 ^, manifest the bilinear character expected as 
their values scale almost perfectly with the coupling constant k = mn. An example is given in 
figure 16.281 where we have carried out four simulations with different values of the perturbative 
parameters. In the upper panel, we have maintained fixed the amplitude of the radial pulsations 
^(bO) _ o.Ol for the H2 overtone, and we have changed the rotation period of the differential 
rotation by an order of magnitude. In the lower panel, the rotation rate is fixed at T c = 100 ms 
and the radial initial amplitude is modified by an order of magnitude. In both cases, the axial 
master function vIa 1 ' 1 ) scales by a factor of ten. Therefore, the results obtained in this section can 
be used to extrapolate the values of the non-linear perturbations for different values of the initial 
perturbative parameters associated with the radial and non-radial perturbations. Of course, this 
procedure is correct as long as the strength of the parameters is within the range of validity of the 
perturbative approach. 

The simulations for the non-linear perturbations are stable for very long time evolutions and 
have a first rate order of convergence, as expected from the implementation of the first order 
accurate Up-wind numerical scheme. 

6.5.4 Effects of radial pulsations on the scattering of a gravitational wave 

In this section, we summarise our study of the scattering of an axial gravitational waves on a 
radially pulsating spherically symmetric star. The aim is to understand whether the outcoming 
signal contains some signature of the radially oscillating dynamics of the star. At first order in e, 
the results are well known and have been re-proposed and briefly discussed in section lrj.4.31 The 
scattered gravitational wave displays an excitation of the first w-mode. and the associated strongly 
damped ringing phase. 

The initial values for the first order axial perturbations are set up for the harmonic index I = 2 
with the condition A6.87I ). The velocity perturbation (3^^ then vanishes and the axial master 
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Figure 6.28: Axial master function ^n 1 ' 1 ), in km, for four different initial values of the first order pertur- 
bative parameters. The upper panel displays the wave forms obtained for radial perturbations excited with 
a H2 overtone of amplitude A^ 1,0 ^ = 0.01 and for axial perturbations with rotation periods T c = 10 ms 
(solid line) and T c — 100 ms (dashed line). In the lower panel the rotation period is fixed at T c = 10 ms 
and the amplitude of the radial pulsations is changed by an order of mag nitude, = 0.01 (solid line) 

and = 0.001 (dashed line). 

function v^ ' 1 ) is given by the Gaussian pulse (16.1051) centered at ro = 20 km with amplitude 
^(O- 1 ) = 0.1 km and width parameter q = 1.25 km~ 2 . The radially oscillating phase is excited 
with the radial eigenf unctions as we have already done in the previous section l6~.5 .31 However, in 
order to have a stronger coupling we have increased the initial amplitude of -y^' ) by an order of 
magnitude, i.e. A^ = 0.01. 

In figures ( I6.29I >- (I6.32I) . we have plotted the first order gravitational master function vl/t ' 1 ) 
and the second order correction due to the coupling vl/v 1 ' 1 ) on a logarithmic scale. The second 
order signal provides small corrections, less than 2 percent when the radial pulsations are excited 
by the F-modes, and less than 0.1 percent for the higher radial overtones. These corrections do 
not modify the properties of the wave forms and spectra of the linear analyses. Furthermore, we 
can notice that the amplification of the signal found in the differential rotating configuration is 
not present in this case. The fundamental difference between this and the previous case, when 
the axial perturbation was related to differential rotation, is essentially due to the difference in 
the properties of the source terms. For the scattering, the source mainly acts for a relatively 
short period ot time, which is given by the travel time of the linear gravitational wave across the 
star. After the wave has been scattered, the first order signal still present in the star decreases 
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Figure 6.29: Scattering of an axial gravitational wave on a radially oscillating star which is pulsating 
in the fundamental radial mode. In logarithmic scale, the upper panel displays the wave forms of the 
first order axial function vl/( ,i) (solid line) and the second order axial function vj/(!,i) (dashed line). The 
lower panel shows again with the solid line the wave forms of vf' ' 1 ) and with the dashed line the function 



according to the long time decay power law (section l6.4.3t . As a result at second order a longer 
ringing phase appears which is anyway well below the first order signal. On the other hand, in 
the coupling between radial pulsations and differential rotation the source terms act periodically 
into the star forever, as this model does not contain the back-reaction. As a result, the source has 
more time to couple with the non-radial perturbations. 

Although the non-linear coupling does not change the linear results, we can study whether 
the Ae wave forms and spectra contain some signatures of the radial pulsating dynamics of the 
star. From ft gures l6~29l and l6~30l we can notice that for dynamical times log [t/ (2M)] ~ 1.70 the 
ringing phase presents a small bump, which is due to a small numerical reflection at the stellar 
surface that appears when the Ae gravitational wave is moving out the star. In particular, this 
effect is caused by the discontinity of the source terms, which vanish in the exterior spacetime. In 
addition, the ringing phase of vl/C 1 ' 1 ) seems to contain some interference effects (see figure l6~32l i. 
This suggests the presence of at least two oscillation frequencies. Therefore, we perform an FFT 
of the signal by selecting the part of the wave form after the numerical reflection. In figure 16321 
the spectra show the presence of the I = 2 w-mode and a non-linear harmonics, whose frequency 
is the sum of the radial and w-mode frequencies. It is worth noticing that this non-linear harmonic 
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Figure 6.30: With the same notation of figure l6^29l the wave forms are related to the scattering of the 
gravitational waves on a radially pulsating star where only the first overtone HI is excited. 



exhibits a strong damping and then the characteristic broad shape of the spacetime modes. How- 
ever, as we have already specified these effects are much lower than the linear results and do not 
modify the spectra of the linear analysis (see e.g.figure l6.15l) . 

The code for the analysis of the non-linear scattering of a gravitational wave on a radial 
pulsating star manifests a first order of convergence. This is expected from the discontinuity of 
^l-'P at the stellar surface, due to the presence of the source terms only in the interior spacetime. 
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Figure 6.31: Scattering of the axial gravitational wave on a radially pulsating star. The solid line in both 
panels refers to the wave form of vl^ - 1 ). In the upper panel the dashed and dotted lines denote the wave 
forms of $( 1 > 1 ) for radial pulsations excited by a F-mode and Hl-mode respectively. In the lower panel 
instead, the dashed, dotted and dot-dashed lines describe the wave forms of v^ 1 * 1 ) for radial pulsations 
excited by H2, H3, and H4 overtones respectively. 
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Figure 6.32: Power Spectrum of the wave form vfrt 1 ' 1 ), which is obtained from the scattering of a axial 
gravitational wave on a radially pulsating star. The radial pulsations have been excited by selecting in each 
simulation a particular radial mode, i.e. the F mode and its first three overtones (see labels in the four 
panels). 



Chapter 7 

Conclusions 



With the aim of investigating the dynamics of non-linear oscillations of compact stars and the 
related gravitational radiation, we have presented in this thesis a formalism and a numerical code 
which enable us to study in the time domain the coupling between the radial and non-radial per- 
turbations of perfect fluid non-rotating compact stars. The formalism for the polar perturbations 
has been worked out in a first paper [89], while the formalism and the applications to axial per- 
turbations are presented in |90|. The applications of the polar perturbations will be presented in a 
future work. 

In order to have a well-defined framework to consider the gauge dependence of linear and 
non-linear perturbations, we have found it very convenient to address this topic by using the 
multi- parameter relativistic perturbation theory introduced in fl~9l 11001 . We have then carried 
out an expansion of the metric, the energy-momentum tensor and the Einstein and conservation 
equations in terms of two parameters A and e, where A parameterizes the radial pulsations, e the 
non-radial oscillations and the Ae terms describe their coupling. The spherical symmetry of the 
stellar background allows us to separate in the perturbative tensor fields the time and radial de- 
pendence from the angular variables, by using the standard tensor harmonic basis. Therefore, for 
any harmonic indices (l,m) the linear and the Ae non-linear perturbations have been described 
by a system of perturbative equations which forms a 1+1 problem, where one dimension is given 
by the time coordinate and the other by the radial (area) coordinate r. The tensor harmonic ex- 
pansion separates all the perturbative fields in two independent classes, the axial (odd-parity) and 
polar (even-parity) perturbations, which are defined by their behaviour under a parity transforma- 
tion. The non-linear perturbative equations and the definition of the gauge invariant non-linear 
quantities have been derived by using the 2-parameter perturbative theory in connection with the 
formalism of Gerlach and Sengupta [43 ], and Gundlach and Martin Garcia [48 74 ]. This formal- 
ism describes generic one-parameter non-radial perturbations of a time-dependent and spherical 
symmetric spacetime in terms of a set of gauge invariant perturbations. Our approach consists 
in exploiting the spherical symmetry of the radial perturbations in order to separate the spherical 
and time-dependent spacetime of the GSGM formalism in a static background, which describes 
the equilibrium configuration, and a first order radially pulsating spacetime. We have then fixed 
the gauge of the radial perturbations and have i) used the GSGM gauge-invariant non-radial e 
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variables on the static background, and ii) defined new second order Ae variables, describing 
the non- linear coupling of the radial and non-radial linear perturbations, which are also gauge- 
invariant for general Ae gauge transformations with the radial gauge fixed. We have then derived 
the evolution and constraint equations for the non-linear Ae perturbations both for the polar and 
axial sectors. As expected, in the interior the Ae variables satisfy inhomogeneous linear equations 
where the homogeneous part is governed by the same linear operator acting on the first order e 
non-radial perturbations, while the source terms are quadratic and made of products of A and e 
terms. In the exterior the sources vanish and there is no direct coupling, and the whole dynamics 
is described by the Ae order Zerilli and Regge-Wheeler functions, respectively for the polar and 
axial perturbations. Thus the effect of the coupling is transmitted from the interior to the exterior 
through the junction conditions at the surface of the star. 

We have given a brief discussion of the boundary conditions, focusing on those at the stellar sur- 
face. In order to avoid negative values of the mass energy density near the surface due to the 
Eulerian description of the radial and non-radial polar perturbations, we have adopted an approx- 
imation already used in literature HI 01 1 11021 ISTIl . i.e. we have removed the outer layers of the 
neutron star, which corresponds to neglecting less than one percent of the total gravitational mass 
of the star. This approximation leads to a description of the second order Ae gravitational radiation 
which is accurate to better than five per cent. 

The perturbations at first order have been studied with the GSGM quantities. However, in 
some cases we have redefined these variables or changed them for computational purposes. The 
radial perturbations have been described by four perturbative fields, two metric and two fluid, 
which obey three first order in time evolution equations and two constraints, as there is a sin- 
gle radial degree of freedom. This allows us to set up either a hyperbolic-elliptic formulation 
(HEF) or a purely hyperbolic formulation (PHF). The numerical simulations have shown a good 
accuracy in both cases. The perturbative equations for the polar e non-radial perturbations and 
for the polar Ae coupling rely on the hyperbolic-elliptic system of three equations already used 
in ll8Tll . where the two hyperbolic equations describe respectively the gravitational wave and the 
sound wave, while the elliptic equation is the Hamiltonian constraint. We have found that this 
system of equations is more suitable for numerical integration than others available in literature, 
as the Hamiltonian constraint is used for updating at any time step one of the unknowns of the 
problem. Therefore, the errors associated with the violation of this constraint are automatically 
corrected. However, the sound wave equation is only known for first order polar non-radial per- 
turbations. Thus, in order to determine it at the Ae perturbative order, we have first obtained this 
wave equation for generic non-radial perturbations on a barotropic, time dependent and spheri- 
cally symmetric background. Then, we have used the 2-parameter expansion and determined the 
source terms relative to the Ae perturbations. The axial non-radial perturbations for the linear e 
and non-linear Ae orders have been studied with a system of two equations, the axial master wave 
equation and a conservation equation. The former describes the only gauge invariant metric vari- 
able of the axial sector, while the latter the axial velocity perturbation. At first perturbative order, 
the stationary character of the axial velocity allowed us to study separately the dynamical degree 
of freedom of the spacetime and its stationary part. Linear axial metric perturbations describe 
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the dragging of the inertial frame, the linear axial velocity perturbation represents a stationary 
differentially rotation. 

Since the linear axial velocity induces a stationary differential rotation, the related metric 
solution describes the dragging of the inertial frame. In the axial sector, we have redefined the 
axial velocity in order to have a vanishing quantity on the stellar surface. This definition leads to 
simpler perturbative equations and a better behaviour of the Ae second order source term near the 
surface of the star. 

We have then presented the numerical code for investigating in the time domain the evolution 
of the axial Ae perturbations. This code is based on finite differencing methods and on standard 
explicit numerical schemes for integrating partial and ordinary differential equations. The struc- 
ture of the code reflects the hierarchy of the perturbative frameworks. Starting from the TOV 
solutions and two independent initial conditions for the radial and non-radial perturbations, the 
code evolves at any time step the linear radial A and axial non-radial e perturbations, thus up- 
dates the source terms of the Ae perturbative equations, and eventually simulates the evolution 
of the non-linear perturbations. In order to have more accuracy in the description of the radial 
pulsations near the surface, we have increased the resolution in this region by adopting the fluid 
tortoise coordinate, which is necessary only for the radial perturbative equations. Therefore, we 
have introduced in the numerical code two meshes for describing the radial and non-radial pertur- 
bations separately at first order. In order to have source terms calculated at the same grid point, we 
have also introduced an interpolation that at any time step connects the simulations of the radial 
pulsations between the two meshes. 

In this thesis the static equilibrium configuration is determined by solving the TOV equations 
for a polytropic equation of state with adiabatic index T = 2 and k = 100 km 2 , and for a central 
mass energy density p c = 3 x 10 15 g cm~ 3 . The star then has M = 1.26M© and R = 8.862 km. 
The structure of the radial and axial perturbative equations enables us to set up two independent 
initial configurations: i) scattering of an axial gravitational wave on a radially pulsating star, ii) a 
first order differentially rotating and radially pulsating star. 

The initial configuration for the radial pulsations has been excited by selecting specific radial 
eigenmodes. We have chosen an origin of time such that the radial eigenmodes are described 
only by the eigenfunctions associated with the radial velocity perturbation j0->°). To this end we 
have first determined the wave equation for this variable, then we have managed it in order to set 
up a Sturm-Liouville problem. This eigenvalue problem has been solved with a numerical code 
based on the shooting method, where the eigenfrequencies of radial modes are accurate to better 
than 0.2 percent with respect to published values. We have performed some tests for the part of 
the code dedicated to the radial pulsations. The simulations for any initial radial mode satisfy 
with high accuracy the Hamiltonian constraint and are stable for very long evolutions. The radial 
spectrum, which has been determined with a FFT of the time evolution, reproduce the published 
results with an accuracy to better than 0.2 percent. 

For the first order axial perturbations we have excited the dynamical degree of freedom of 
the spacetime with a standard gravitational wave scattering on the star, where the impinging axial 
gravitational wave is described by a Gaussian pulse. The tests carried out on the numerical code 
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show that we are able to reproduce the wave forms and the spectra known in literature for the 
harmonic indices I = 2, 3. The axial differential rotation instead has been described by expanding 
in vector harmonics the relativistic j-constant rotation law, then taking the first component which 
is related to the gravitational wave emission, that is I = 3. We can specify two parameters in the 
initial profile for the axial velocity perturbations, i.e. the differential parameter A and the angular 
velocity at the rotation axis il c . The value for A has been chosen in order to have a smooth profile 
and a relatively high degree of differential rotation, as for high A the rotation tends to be uniform 
and then the I = 3 component vanishes. We have chosen an angular velocity that corresponds 
to a 10 ms rotation period at the axis of the star. For other values of the angular velocity, the 
linearity of the perturbative equations allows us to get the respective gravitational signal with a 
simple re-scaling. The relativistic j-constant rotation law contains a term which is related to the 
dragging of the inertial frame. In order to simplify the expression of our initial condition we have 
estimated the amount of correction associated with the I = 3 component of this relativistic term. 
We have then neglected it, thereby introducing an error of less than ten percent. 

For the first initial configuration, i.e. the scattering of an axial gravitational wave on a radially 
pulsating star, the second order signal provides small corrections, less than 2 percent when the 
radial pulsations are excited by the F-modes, and less than 0.1 percent for the higher radial over- 
tones. These corrections do not modify the properties of the wave forms and spectra of the linear 
analyses. The second configuration, where the linear axial perturbations describe a stationary 
differential rotation and the associated frame dragging, produces an new interesting gravitational 
signal. The wave forms have these properties: i) an excitation of the first u>-mode at the early 
stage of the evolution ii) a periodic signal which is driven by the radial pulsations through the 
source terms. The spectra confirm this picture by showing that the radial normal modes are mir- 
rored in the gravitational signal at non-linear perturbative order. However, the excitation of the 
u>-mode at the early stages of the numerical simulations is an unphysical response of the system 
to the initial violation of the axial constraint equations for the coupling perturbations. Moreover, 
a resonance effect is present when the frequencies of the radial pulsations are close to the first 
u>-mode. For the stellar model considered in this thesis the amplitude of the gravitational wave 
signal related to the fourth radial overtone is about three orders of magnitude higher than that 
associated with the fundamental mode. We have also roughly estimated the damping times of the 
radial pulsations due to the non-linear gravitational emission. Their values radically depend on 
the presence of resonances. For a 10 ms rotation period at the axis and 15 km differential pa- 
rameter, the fundamental mode damps after about ten billion oscillation periods, while the fourth 
overtones after ten only. This is not surprising, and shows that the coupling near resonances is a 
very effective mechanism for extracting energy from the radial oscillations. 

It is worthwhile to remark that a possible detection of this gravitational signal could provide 
new information of the stellar parameters, as the second order gravitational spectra reproduce 
those of the radial modes of a non-rotating star, which can be determined easily for a large class 
of equations of state. The numerical code manifests a first order convergence for the simulations 
of the non-linear axial perturbations arising from the coupling between the first order radial and 
axial non-radial perturbations. 
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Future extensions 

The implementation of the numerical code for studying the coupling between the radial pulsations 
and the polar non-radial perturbations is currently under way. The spectrum of the polar non-radial 
perturbations is richer than the axial case due to the presence of the fluid modes, which may have 
frequencies lower than the spacetime modes and a longer gravitational damping. As a result, we 
may expect a more effective coupling with the radial pulsation modes. 

Future extension of this work certainly must consider more realistic models of the star, which 
take into account the effects of rotation, composition gradients, magnetic fields, dissipative ef- 
fects, etc. In particular, it would be interesting in a protoneutron star to compare the damping rate 
due to the gravitational radiation produced by the coupling with the strong damping induced by 
the presence of a high-entropy envelope, which surrounds the newly created neutron star. 

New interesting non-linear effects could come out when the star is rotating, mainly due to the 
different behaviour of the radial and non-radial modes in a rotating configuration. While the radial 
modes are only weakly affected and their spectrum can be considered almost the same as that of 
a radially pulsating non-rotating star (when scaled by the central density), the non-axisymmetric 
modes manifest a splitting similar to that observed in the atomic energy levels due to the Zeeman 
effect. The rotation has the effect of removing the mode degeneracy of the azimuthal quantum 
number, which is present in the non-rotating case. The amount and the details of this frequency 
separation depend on the stellar compactness and rotation rate. We might then expect that for a 
given stellar rotation rate and compactness the non-radial frequencies should cross the sequence 
of the radial modes l331l . In this case the frequencies of these two kinds of modes are compara- 
ble, and possible resonances or instabilities could influence the spectrum and wave profile of the 
gravitational wave radiated. The possibility of identifying new resonances in the high frequency 
gravitational wave spectrum could provide new relations, which can be used to determine the 
stellar parameters through asteroseismology. 



Appendix A 

Gundlach-Garcia Source terms 



In this section we re-write the source terms of the equations of polar non-radial perturbations on 
a time dependent and spherically symmetric background, which have been found by Gundlach 
and M. Garcia [48 ]. These equations are valid for a perfect fluid star with constant entropy along 
each element fluid trajectory. The equations for a barotropic fluid can be easily determined by 
neglecting the terms relative to the entropy s, its perturbation a and the quantity C = ( §f 
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In this section, we report also the equations relative to the radial perturbations of a time dependent 
and spherically symmetric star, which have been determined in the radial gauge by Gundlach and 
M. Garcia [48]. The two fluid evolution equations are given by 
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Appendix B 

Sound wave equation 



Here, we give the complete form of the sound wave equation (I4.36t for non-radial perturbations 
on a barotropic, time dependent and spherically symmetric background. This equation is written 
in terms of the quantities and frame derivatives introduced in the GSGM formalism i2.2l . 
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Appendix C 



Source terms for the Ae polar 
perturbative equations 



The perturbative equations which describe the coupling between the radial and non-radial polar 



perturbations (I4.37b - (l4.42b have long terms that are written in this appendix. For the gravitational 



wave equation (I4.38t the source has the following form: 
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' + 1S r(3-(A,+«,)r)^ 



f_ 7 (i.o) 
r 



-$-2A 



(C8) 



The source of the sound wave equation ( I4.39t is given by ( 14.401 ): 

S H = hHW + 62jff (o/) + 63^.1) + + (o,D + M (p4) + ^5(0,1) 

fc (o,i) _ ^5(0,1)) 1 + 6g ^5(0,1) + fc (o,i)) + bl()7 (o,i) + feu7 (o,i) + 6l2 ^(o,i) 



+ &8 



(C.9) 



and the expression of the coefficients 6j is the following 



bi 
h 
h 



+ 



+ 



47rr 



dp 



-2A 



2 (1 - c 2 ) e 



„-*-A (1,0) 
-2A 



r -» A P + 0, r dcf + 2c 2 p-p \ / 2 (l,0) e *' 
r; r/p p / V 



2tt c 2 dp 



+ (A )P + $ r ) 



e- 2A A P + $ r dc 2 



27rr 



— - C 2 + 1 + l) 7 d.0) e »-A + ^(1.0) I g -2* 



(CIO) 
(C.ll) 

(C.12) 



1 

47rr 



A jP I A, 



1 



+ $ r (2$ r + A P ) 



2(l-c 2 )^°)+2 



2$ r - A r + - ) c 2 - <i> 

r 



1 dc 2 
c 2 dp 

(1,0) 



F (1.0) e -4A 



rp 



3$ r + 1 + (A, r - 4*, - 1 ) c? 
-2A 



. 5 (1,0) I e -2A 



(C.13) 



^ ( 1 - e 2A + *, r A, r r> + ( A P + ^ A ' + ■ * ^ - ^ 



dp 



(A jP + $ P ) (1 + 3c 2 ) 



K[±l)-2 e 2A 



5(1,0) 



+ 



c 2 r 



1 ^ 66 2 + 1+r$r) A J r+^r e -2A 



4-7rr 



+ dc| 
dp 



+ 3 C 1 + 4c 2 + l]id^) + e " 2AA - + ^ dC " : 



4irr 



dp 



--A, P + 2*, r fT^+^°) 
r 



+ - 



(A r + $ r )(l + 3c- 2 )-^ +1) ^ 2A 



(A, P + $, r ) h+^_£j^_(i_^)$. 



^(1,0) I, g -2A 



(C.14) 



7 (i,o) e * 



A P + - c 2 + T 



7 (i,o) e # _ £ 



r 2 7 d,0) e * 
T" 2 V / ,r 



-2$-A 



(C.15) 
(C.16) 
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+ 



bio 



+ 



bi2 = 
bi3 = 



+ 



+ 



87rr 



dp 



i7 (l,0) e -2A + (1 _ -2) ^ ^(1,0) + l H £>°A 



-2A 



(C.17) 



J_ (l- e 2A + 3r* r + rA, r )^-4*, 



,(i,o) 



e 2A -l 



r 

A.r + 



2$ r (4 + rA >r ) - 3A 

2 ' 2A 1 rffr 
1 + 2$ r r + 3c 2 



c 2 + 4r$ 2 r 



5(1,0) 

((3$ r + A r ) r + 1 - e 2A ) 



47rr 2 c 2 dp 

+ 2 2$ r + (A, r - 2$ r - ^ c 2 i^°)-2c^£°>} e - 2 \ 



jy(i.o) 



(C.18) 



2 4 



( r 2 7 (1.0) e *) 



2p 



2 (1,0) * 



-$-2A 



(C19) 



1 



^ < 2rc 2 (1 - c 2 ) (W^e*) ^ + [|j (A, r + *,) ^e- 2A + 2 (*, + 2A r + r g 



-.1 



4(l + rA, r ) + ^(2 + r$ r ))c 2 



2 - - J r$. 

5 $. 



(r 2 7 (1 ' 0) e*) 



+ I -2 (A r + $ r ) (1 + (* r - A r ) r) 4 + - ^ (2$ r r + 7 + 8r A 

,2A 



- 2A, r (rA, r - 1) + 



p $ r r (2$ r r + 5) + 2 e 



2r 



2r 



1 + |)(2+ <!>.,,■) )r 



,-2A 



4 - |J (A r + * >r ) r + 4) * >r - ^ (A, r + d>, r ) 2 * r ^| 



-$-2A 



C s 

2r 2 
1 



(A, r + *, r ) (1 + c 2 ) + t (rV^e*) - (1 + 



-*-2A 



(C.20) 
(C.21) 



2 1 « 



(1+ 



-2A 



r ^ (1 ' 0)e$ ) r + l^T (V + * r) *r [(A,r + *r) ( 



r 2 7 (l,0) e * 



( r 2 7 (1.0) e * 



dp 



(2A r + *, r ) + ( - r + 2A r + 3* r + V - ^ - ^ 



c + 2 



2p ' ' 



( r 2 7 (1.0) e *) 



4 A r + $ r 

- + — ^ ^ 

r 2 



6 



p 



(A >r + * r ) A >r - $ r + 



1 

2 
1 

4 r 2 
2 



6A r + 7$ . 



p 



P 



2A 



r 2 7 (l,0) e * 



+ 



-*-2A 



A. 



2r 



A. r + 



(C.22) 



614 = — 



p J p r 



(r 2 7 (1 ' 0) e*) 



- (A^ + ^^rV^e* 



-2A 
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The third equation of the system is the Hamiltonian constraint, whose source term has the 
following form ( 14.421 ): 



Snamii = * (kW> - S$») + c 2 k^ + C3k fV + C45 (o,D + C5k (o,D + CqH 



+ c 7 ^jf l) + c 8 V iU,ij + c 9 7 

where the coefficients Cj are given by 

r 5(1,0) 



(0,1) 



,(o,i) 



Cl 
C2 

C3 
C4 

C5 
C6 

C7 
C8 
CQ 



3^(1,0) + A,r + $,r g (i, )^ j 

(A ir + $ r)e A-* 7 (i,o) ; 

5(1,0) + 2 
2 A r + $ 



2 A r + $ . 



1 + 



1 e~ 2A A r + $ r dc 2 



47rr 



dp 



r 

l r 



(2 - 4A, r r + /(/ + !) e 2A ) 7 (1 ' 0) + 2r 7 



-^(A, r + $ r ) 7 ( 1 .o). 



,(i.o) 



(0,1) 

(C.24) 

(C.25) 
(C.26) 

(C.27) 
(C.28) 

(C.29) 
(C.30) 

(C.31) 
(C.32) 
(C.33) 



(C.23) 



Appendix D 



Source terms for the Ae axial 
perturbative equations 



The equations which describe the coupling between the radial and non-radial axial perturbations 
are given by the two equations (I4.62t - d4.63l l. Their source terms and have the following 
form: 



= 2 (r^ 1 '°)-r ? ( 1 ' )) e- 2A *(o,i) + | 4 ^ ( ^ + p )r 1^^(1,0) 



+ 



, 9 4M 
47r (p - p) r A - H 



5 (i,0) + 2 



'a f- —\ 2M 
Air (p- p)r 2~ 



,,(1,0) I $(0,1) 



7 (1.0) +87r (p + p) re A+* 7 (l,0) 



(0,1) 



4vr (p + p)L^L H W 



+ 



' , , 1(1 + 1) + 3 , 12M 
4tt (p - p) r - -i '- + — =- 



5(1,0) 



+ 2 



' + 6M 

4vr (p - p + - - — 



^(1,0)1 $(0,1) _ 8 ^ r ( 4r? (l,o) _ 3r5 (i,o)) e-A^o.i) 



+ 1 24vr (4vr pr 3 + M) e A S (1 ' 0) - 32vr ^4vrp r 2 + — e A r) 



A„(1,0) 



+ 167rre- A ^°) 



,1) 



(D.l) 



^ = - 7 a.o) e -A^(o,D + 



M 

An pr J 



,2A 



t \ 7 (i.o) 



~(i.o) 



-0(0,1) 



(D.2) 



Appendix E 

Tensor harmonics 



In this section we write the even parity (polar) and the odd parity (axial) tensor harmonics, which 
have been defined in a covariant way in section l2~2.2l The polar scalar spherical harmonics, which 
satisfy the condition A2.49I ). are given by the following explicit expression: 

Y lm (0, <j>) = e im \l 2 \ +1 f~ n f ] Pi m (cos 6) , (E.l) 
y Air [l + m)\ 

where P\ m are the associated Legendre functions with the harmonic indices (l,m), which are 
given by 

P lm (*) = {-IT (1 - x 2 ) - J^P (x) = L^- (1 - X 2 ) - |^ (x 2 - 1) ', (E.2) 

where x € [—1,1], and in the last equality we have used the "Rodrigues representation" for the 
Legendre function Pi(x): 

flM-yflSr^- 1 )'- <E3) 

It is worth noticing that for axisymmetric stellar perturbations m = 0, the associated Legendre 
polynomial Pi m reduces to the Legendre polynomial p. 

The polar and axial vector harmonic bases are then given by 

ylm = ylra = (yin yl^ ^ (R4) 
S lrn = e bylm = (_^_ylm }Shld yln i y (R5) 

The tensor harmonics for the polar sector are given by 

vim ~V~lm„, rylm vim , -Q \^lrn„ J /ri r\ 

Y ab = Y lab: ^ab = Y :ab H 3 7a& ' ^ ^ 
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where the explicit expression of the second order covariant derivative of the spherical scalar har- 
moncis is given by the following expressions: 



vim 



Y% Y%- cot 6Y% < 

Y l Z - cot 9 Y l T Y l Z + sin 9 cos 9 Y lr ! 1 



V 

The axial sector has the following tensor harmonics: 

aim — olm _i_ aim 
a ab = a a:b + ^b-.a ) 

where S l J£ assumes the following form: 



(E.7) 



(E.8) 



/ -- L 



~ilm 
V b 



Y l Z. - cot 9 Y 



Irn 



1 vim 



Y l Z, - cos 9 y' fl m sin 9 ( Y l Z — cot Y 



sin#Y 



Zm 



/Im | 



V 



(E.9) 



/ 



The harmonic tensors form an orthogonal basis of the 2-dimensional sphere S 2 . The orthogonality 
relations for the polar sector are defined as follows: 



vim \ yi'm' 



J dn ( 



gll' gmm' 



IV rmm' 



(E.10) 
(E.ll) 

y dft 7 a S cd ) * Z& m ' = 1(1 + 1) 12 V '<5 mm ' , (E.12) 



/ (2 + 1) 5" <5 



where S u> is the Kronecker delta and the asterisk denotes complex conjugation. For the axial 
sector the orthogonality conditions are given by 



J dfi 7 a6 (S l a m ) * Sl m ' = 1(1 + 1) 5 u 'S mrn ' , (E.13) 
J dQj a \ cd (S l { ™ c) y S( b % = 21(1 + 1) (l 2 + I - 1) b ll 'b mm ' , (E.14) 

where we have denoted (a : b) = a : b + b : a. 



Appendix F 

Finite difference approximations 



u j+ i = Uj + Axu'j + -^u" + 0(Ax 2 ) , (F.2) 



In "finite differencing method" the derivative operators are approximated by finite difference for- 
mulae. The finite approximations are based on an appropriate linear combination of the Taylor 
expansions of the function of interest. In this section we illustrate this procedure only for the first 
derivative u' = u )X of a one-dimensional scalar function u = u(x). For higher derivative orders 
and variable numbers the method is similar. Let u = u(x) be discretized on a one-dimensional 
mesh of dimension J and increment Ax, 

uj = u(xj) for j = 1, .., J . (F.l) 

The Taylor expansion of Uj + \ and Uj-i around the grid point Xj reads: 

Ax 

i 

2 

Ax 

= Uj - Axu'j + —u'j + 0(Ax 2 ) . (F.3) 
Now, when we take their difference and isolate the term u'j we obtain: 

Don = U]+1 ~^ 1 + 0(Ax 2 ) , (F.4) 

where Dqu denotes the centered finite difference approximation of the derivative u'j , which shows 
a second order accuracy 0(Ax 2 ). For grid points near the surface these centered formulae are 
not able to approximate derivatives by using only internal points of the mesh. For this aim it is 
more appropriate to use one-sided approximations of u'. First order accurate approximation of u' 
is given by the following expressions: 

D+u = w-^ + 0(Ax) (R5) 
Ax 

D u = U 3~ U 3-1 + 0(Ax) (R6) 

Ax 

which can be determined by the Taylor expansions dF2t and dF3l >. 
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Second order accuracy is instead reached with the following finite difference formulae: 

D + u = _ 3u,--4y+u i+2 + 2 

= —2 3 —± J — + 0(Ax 2 ). (F.8) 

The finite expressions of the second derivative u" can be determined on the same line as the 
first order. The centered finite approximation is given by the following expression: 

D ^ u j+1 -2n,+ Uj _ 1 +0{Ax2h (R9) 
while the one-sided formulae for u" are the following: 



D+u = ^ J ^ Ax2 ^ ^ +0(Ax 2 ), (F.10) 
„o 2m — 5u,_i + 4uo_9 — it,-s „ , . 



Appendix G 

Numerical methods 



The numerical algorithms used in the simulations are all present in the reference "Numerical 
Recipes" l9T1l . Here we briefly illustrate the McCormack scheme that is not mentioned in this 
reference. Furthermore, we also write the methods used for determing the convergence rate and 
the norms of the perturbative solutions. 



G.l McCormack algorithm 

It is an explicit second order differencing algorithm and a 2-level method. The scheme consists of 
two computational steps, namely a predictor and a corrector step. For a given partial differential 
equation of a function / = f(t, x) that we indicate as follows: 

f,t = G(f,f, x ,t,x) , (G.l) 



fn+l =AtG n (G _ 2) 



the value of the function / on the new time-slice is updated as follows 
i) Estimator step 

where the source term Gj is evaluated at £j_i/2 by using the values fj l __ l and /j 1 . 
i) Corrector step 

where now the term G]^ is determined by the preliminary values /j 1 and /j+i . 



G.2 Convergence test 

The exact solution of a given analytical equation has to be more and more accurately approxi- 
mated by numerical solutions determined with an increasing resolution of the numerical mesh. In 
the limit of an infinite dimension of the grid the numerical solution must tend to the exact solution. 
In addition, the rate of the convergence has to be consistent with the degree of approximation of 
the numerical method used. 
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Let f ex be an exact solution of a given equation, and f c and f m the numerical solutions 
determined respectively on a coarse grid of dimension J c and a medium grid of J m = 2 J c . The 
deviation of these solutions from the analytical one can be then written as: 



0(r 



^T+l^ 



A/™ EE 



^g(/r-/f) 2 =^)W- 



(G.4) 



(G.5) 



where the differences in these two expressions are both evaluated at the points of the coarse mesh. 
The letter a denotes the accuracy order of the numerical solution and Eq is the unknown error 
term. The ratio of the previous two expressions leads to the following expression: 



AT 

Af m 

thus the convergence factor a is given by 



= 2< T + 0(r <T+1 ), 



(G.6) 



a = 



log [A/7 A/" 



log 2 



(G.7) 



When the exact solution is unknown the convergence test requires three numerical solutions 
determined on three different grids whose resolution is in the following proportion 1:2:4. There- 
fore, in the previous expressions we can replace the exact solution f ex with the numerical solution 
f* obtained on the fine grid of dimension Jf = 2J m = 4 J c . It is important to remark that in this 
second case the convergence test informs us of the correct scaling of the numerical error but it 
does not imply that the numerical solution is going to converge to the true solution. 



G.3 Numerical stability and dissipation 

The numerical stability and dissipation of the simulations can be monitored by the constancy of 
the norms, which can be determined for numerical solutions during their time evolution. Let /" 
be the finite approximation of a quantity /, which has been determined on a one-dimensional 
grid of dimension J at the n time slice. The L2 norm can then be calculated at any time step as 
follows: 

H/ll2 = 7E(/™) 2 > (G.8) 
J j=i 

where we have introduced the division by J for having an expression averaged on the number of 
grid points. 
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